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COMPUTER  MODELS  FOR  TWO-DIMENSIONAL 
STEADY-STATE  HEAT  CONDUCTION 


M.R.  Albert  and  G.  Phetteplace 


INTRODUCTION 

Most  major  Army  facilities  are  heated  by  central  heat  distribution  systems.  Heat  losses  from 
these  distribution  systems  are  of  interest,  as  arc  temperatures  of  soil  surrounding  the  distribution 
pipes,  for  several  reasons.  If  we  are  to  design  efficient  heat  distribution  systems  we  must  be  able 
to  accurately  assess  the  heat  losses  in  order  to  determine  the  optimum  trade-off  between  the  cost 
of  insulation  and  the  continuing  cost  of  heat  loss.  In  instances  where  portions  of  a  system  have 
failed  and  require  replacement,  again  we  need  to  accurately  assess  potential  heat  losses  to  determine 
the  optimum  insulation  thickness  for  the  areas  that  must  be  replaced.  In  areas  of  seasonal  frost 
and  permafrost,  the  temperature  distribution  in  the  soil  around  buried  heat  carrying  pipes  is  of 
primary  importance.  Significant  thawing  of  permafrost  can  cause  loss  of  support  and  subsequent 
pipe  settlement  and  possible  pipe  breakage.  In  areas  of  seasonal  frost,  piping  systems  should  be 
buried  deep  enough  to  prevent  freezeup  in  the  event  of  a  system  outage. 

Several  methods  exist  for  analyzing  heat  transfer  problems  of  this  type.  Some  of  the  methods 
which  have  been  applied  to  conduction  heat  transfer  include:  1)  analytical  methods,  2)  approxi¬ 
mate  methods,  3)  empirical  methods,  and  4)  numerical  methods. 

Analytical  methods  are  useful  for  a  limited  class  of  problems  for  which  closed  form  analytical 
solutions  to  the  heat  conduction  equation  have  been  found.  They  arc  limited  because  only  the 
simplest  geometries  may  be  treated  with  corresponding  ideal  boundary  conditions.  One  such  solu¬ 
tion  exists  for  a  buried  uninsulated  single  pipe  in  a  semi-infinite  medium  with  uniform  and  constant 
material  properties  and  boundary  conditions;  this  will  be  discussed  later.  Because  of  the  very  limited 
class  of  problems  for  which  closed  form  analytical  solutions  have  been  found,  they  have  little  appli¬ 
cation  to  real  world  engineering  problems.  However,  they  can  be  used  to  find  approximate  solu¬ 
tions  to  actual  problems. 

Approximate  methods  are  actually  an  extension  of  analytical  methods.  In  some  instances  an 
exact  solution  to  the  heat  conduction  equation  cannot  be  found,  yet  by  making  some  reasonable 
assumptions  an  approximate  solution  may  become  possible.  Approximate  solutions  arc  limited  to 
certain  problems,  like  analytical  solutions,  except  to  a  lesser  extent.  They  can  provide  estimates 
suitable  for  applications  not  requiring  precise  results. 

Empirical  methods  are  based  on  experimental  results  instead  of  on  the  solution  of  the  governing 
equation.  With  knowledge  of  the  form  of  the  governing  equation  and  solution,  experimental  data 
can  be  used  to  find  empirical  equations  which  approximate  the  physical  process.  A  typical  example 


of  this  approach  is  the  development  of  conduction  shape  factors  by  electrical  analog  experiments. 
Again,  empirical  methods  are  limited  to  cases  where  empirical  equations  have  been  found.  These 
equations  are  usually  only  applicable  to  very  specific  problems  and  often  cannot  be  accurately 
extrapolated  to  describe  similar  problems. 

Numerical  methods  have  found  many  applications  in  heat  transfer.  In  short,  a  numerical  method 
uses  an  approximation  of  the  governing  differential  equation,  or  its  solution,  at  a  number  of  dis¬ 
crete  intervals  of  a  region  to  find  an  approximate  solution  to  the  heat  conduction  equation  over 
that  discretized  region.  The  enormous  advantage  of  numerical  methods  is  their  ability  to  model 
nearly  any  problem,  including  ones  without  ideal  boundary  conditions.  This  is  an  advantage  no 
other  method  can  claim.  The  disadvantage  of  numberical  methods  lies  in  their  complexity.  Cri¬ 
teria  governing  the  stability,  convergence,  consistency,  and  accuracy  of  a  numerical  method  must 
be  established  before  the  method  may  be  considered  valid,  and  for  anything  but  the  simplest  prob¬ 
lems  a  computer  is  required  for  solutions.  Fortunately,  however,  the  computer  program  can  be 
written  to  handle  a  broad  class  of  problems  and  the  user  need  only  be  familiar  with  certain  aspects 
of  the  problem  to  accurately  use  the  program  to  find  the  results. 

The  purpose  of  this  report  is  to  outline  the  finite  difference  and  finite  element  numerical  method 
as  they  apply  to  steady-state  heat  conduction,  to  compare  the  results  of  each  method  to  analytical 
results  and  also  to  the  results  of  the  other  method,  and  to  demonstrate  the  use  ot  each  method’s 
computer  program  to  the  potential  user. 


THE  FINITE  DIFFERENCE  METHOD 


Basic  ideas  of  finite  differences 

The  finite  difference  method  is  a  numerical  method  used  to  approximate  the  solution  of  a  par¬ 
tial  differential  equation.  Common  partial  differential  equations  for  which  finite  differences  have 
been  employed  include  the  wave  equation,  the  vorticity  equation,  Poisson’s  equation,  and  the  heat 
conduction  equation, 


Often  boundary  conditions  or  nonlincaritics  in  a  problem  involving  differential  equations  com¬ 
plicate  the  problem  so  that  analytical  solutions  are  unavailable  and  numerical  solutions  must  be 
used.  The  finite  difference  method  approximates  each  partial  derivative  in  the  differential  equation 
by  an  algebraic  finite  difference  representation. 

Finite  difference  representations  may  be  derived  either  from  a  Taylor  series  expansion  about  a 
point  or  from  physical  considerations,  such  as  an  energy  balance  in  the  case  of  heat  conduction. 

Let  us  briefly  review  the  Taylor  series  procedure.  Consider  the  forward  Taylor  expansion  of  a 
function  f(x)  about  x: 


f(x  +  Ax)  =  f(x)  +  Axf’(x)  +  tm(x)  +  ...  (I) 

l  o 

where  the  primes  denote  differentiation  with  respect  to  x.  This  may  be  solved  for  /'*(*)  by  trun¬ 
cating  the  higher -order  derivatives  to  obtain 

.0|M,  (2) 

This  is  the  forward  finite  difference  expression  for  the  first  derivative  of  the  function  1  (x).  The 
term  O(Ax)  indicates  that  the  error  due  to  truncation  of  the  Taylor  scries  is  ol  the  order  Ax. 


The  representation  is  commonly  used,  for  example,  as  an  approximation  of  bTjbt  in  the  heat  con¬ 
duction  equation,  where  we  let 


bT  =  T(t  +  At)  -  m  (3) 

bt  At 

where 

T(t)  -  temperature  at  time  t 
At  -  time  step. 

Similarly,  we  may  find  the  backward  finite  difference  expression  for  the  first  derivative 
f(x  -  Ax)  =  f  (x)  -  (Ax)f’  (x)  +  f"  (x)  -  ^  f"  (x)  +  ...  (4) 


where 

f'(x)  =  ~  — )■  +  O(Ax) . 

To  obtain  a  difference  representation  for  the  second  derivative,  eq  1  and  4  arc  added 

f{x  +  Ax)  +  f(x  -  Ax)  =  2f(x)  +  (Ax)2f”(x)  +  0[(Ax)2] 

ftx\-  f{x  +  Ax)  -  2f(x)  +  f(x  -  Ax) 

(Ax)2 

Equation  5  is  the  central  difference  representation  for  the  second  derivative.  Note  that  it  is  accu¬ 
rate  to  the  order  of  (Ax)2.  It  is  evident  that  the  accuracy  of  a  solution  found  using  finite  differ¬ 
ences  will  be  dependent  upon  our  taking  Ax  sufficiently  small. 

When  using  the  finite  difference  method,  we  set  up  a  network  of  grid  points,  called  "nodes,”  in 
the  region  of  the  spatial  independent  variables  in  which  we  are  interested.  The  boundaries  of  the 
grid  should  coincide  with  those  of  the  problem.  The  partial  derivatives  of  the  original  partial  dif¬ 
ferential  equation  are  replaced  by  their  corresponding  difference  representations,  and  an  equation 
is  set  up  for  each  point  in  the  grid.  The  resulting  set  of  algebraic  equations  may  then  be  solved. 

Finite  difference  computer  program 

A  finite  difference  computer  program  was  written  to  model  steady-state  two-dimensional  heat 
conduction.  The  program  has  been  set  up  in  a  general  form  to  allow  many  different  conductivities 
in  the  region  of  interest  and  to  permit  constant  flux,  convective,  constant  temperature,  and  semi- 
infinite  boundary  conditions.  The  general  form  of  the  program  allows  new  problems  to  be  set  up 
very  easily. 

Application  of  finite  differences  to  steady-state  heat  conduction 

For  two-dimensional  steady-state  heat  conduction,  heat  transfer  obeys  the  elliptic  partial  differ¬ 
ential  equation 

where 


T  -  temperature 
k  -  thermal  conductivity 
x,y  =  spatial  variables. 
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If  the  conductivity  k  were  to  be  constant  over  the  region,  the  finite  difference  form  for 
d/dx(kdTlbx)  =  fcd2  779x2  would  be  immediately  obtainable  from  eq  5, 

k  &T  =  r(x  +  Ax)-2r(x)  +  7~(x-Ax) 
dx 2  (Ax)2 

However,  the  finite  difference  computer  program  was  to  be  written  so  that  conductivity  could  be 
a  function  of  location.  Then  the  resultant  conductivity  must  be  used  in  regions  of  varying  conduc¬ 
tivity.  The  “resultant  conductivity”  is  easiest  to  understand  by  recalling  that  conductivity  is  the 
reciprocal  of  resistance.  The  net  resistance  between  any  node  x,y  and  node  x  +  Ax,  y  is 


R  = 


Rx  +  R 


x+4x 


°x  °x+Ax 

k  +  k 

"x  "x+Ax 


where,  for  example,  dx  is  the  fraction  of  the  distance  along  the  heat  flow  path  which  is  associates 
with  conductivity  kx.  In  this  program,  space  increments  Ax  and  Ay  arc  assumed  to  be  uniform, 
then  dK  =  tfx+Ax  =  '/a.  The  resultant  conductivity  is  the  reciprocal  of  resistance  R: 

1  -  2 

*x  *x+ix 


Thus  the  finite  difference  formulation  for  each  term  in  eq  6  is  as  follows: 
L  L  in -/  2  \  T(x  +.Ax) -  7~(x)  /  2  \ 

3x  V  ax/  (  _j _ ^  (Ax)2  [_J_  +  _L 

\^x+Ax  *x/  \^x-Ax  kx! 


T(x  -  Ax)  -  T(x) 
(Ax)2 


(7) 


'■y+Ay 


\  T(y  +  Ay)  -  T(y) 

(  2  > 

7~(y- Ay)  -  T(y) 

j  (M2  ( 

1  +  » 
(^y-Ay  ^y/ 

)  (Ay)2 

(8) 


Combining  eq  7  and  8,  we  arrive  at  the  finite  difference  formulation  for  eq  6 


(9) 


Here  it  is  convenient  to  change  the  notation  to  a  more  commonly  used  form.  Let  T-t  j  represent 
the  temperature  of  the  node  under  consideration,  where  /  is  increased  in  the  vertical  direction  and 
/  in  the  horizontal  direction.  The  grid  set  up  will  appear  as  shown  in  Figure  1  (the  length  and  width 
of  the  grid  are  variable,  and  may  be  established  when  data  are  entered  into  the  computer  program). 

Each  node  represents  the  area  (enclosed  in  dotted  lines  in  Fig.  1)  around  it.  The  most  conven¬ 
ient  way  to  set  up  the  computer  program  to  enable  it  to  handle  many  different  cases  is  to  use  a 
square  grid,  that  is,  let  Ax  -  Ay  =  As.  Now  we  multiply  each  term  in  eq  9  by  (As)2,  change  to  the 
new  notation,  and  rearrange  the  terms  to  arrive  at  the  usual  finite  difference  equation  for  a  node 
IJ  not  on  a  boundary, 
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Note  that  each  node  represents  a  region  of  space,  and  will  therefore  have  a  specified  thermal  con¬ 
ductivity.  The  above  equation  allows  for  a  different  conductivity  for  each  node  of  the  grid. 

Boundary  conditions 

The  finite  difference  equations  used  in  heat  transfer  may  be  derived  either  by  replacement  of 
a  partial  derivative  by  its  difference  representation,  as  above,  or  by  calculating  an  energy  balance 
between  a  node  and  each  of  its  surrounding  nodes.  The  boundary  conditions  in  the  following 
sections  will  be  derived  using  energy  balance  considerations. 

Constant  flux  boundary 

For  a  node  on  a  constant  flux  boundary,  the  condition  dT/dx  =  C  exists  at  the  boundary.  Con¬ 
sider  a  node  on  the  right-hand  boundary.  The  control  volume  is  that  which  is  associated  with  the 
node,  as  illustrated  in  Figure  2.  Examine  the  heat  flow  between  node  /',/  and  its  surrounding  nodes. 
Between  node  /,/  and  node  i-\,j  we  have,  for  unit  depth, 
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i - r~ 

•  I  •  1 

i 

-I - 1 - h  - 

1  .  1  .  '  ■ 

I  i  |  T*.i 

■j - 1  —  +  - 

*  !  *  !w 

A* 


>Ay 


Figure  2.  Node  on  boundary  of  finite  difference  grid. 


where 


1  +  _L 

Ax 


resultant  thermal  conductivity  between  nodes  i- 1 ,/  and  i,j 


-  surface  area  perpendicular  to  the  direction  of  heal  flow,  for  unit  depth 


^i-t  ,i  ”  ^i,i 


Ay  -  distance  between  the  nodes 

=  temperature  difference  between  nodes. 


Similarly,  for  the  heat  flow  from  node  /'+1 into  node  /,/  we  have 

2 


-L_  +  _!_ 

.  k.  .. 


2^ 


i+  t ,  j  "i.j 
and  the  flow  from  node  i,  /'-I  into  node  i,j  is 

Av 


ten) 

Ki-1  *u  / 


Ax 


(^i-i-^i)- 


The  heat  flow  crossing  the  boundary  into  node  ij  is  given  by 
f?4  ”  0  '  As 

where  tp  is  the  heat  flux  crossing  the  boundary  per  unit  area.  At  steady  state  the  sum  of  the  heat 
flows  is  zero.  Then  we  have 

Of  +  Q?  +  Q 3  +  Q4  =  0  ■ 

For  Ax  =  Ay  =  As, 

(p=T X)  r»->  *  (nf-r) *  (  J1 ;  Zjj) 
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/ - 5 -  +  - - - + - ! - \  j. . 

(— ! —  +  _L  —1_  +  _L  _J —  +  JL)  *■’ 

\*i-l ,  j  *i,i  * i.  j  *1+1, j  *i,i/ 


=  -<t>&S  . 


This  equation  applies  to  the  right-hand  boundars  ;  the  indices  may  be  appropriately  rearranged  to 
obtain  the  equation  for  other  boundaries.  Note  that,  for  a  boundary  which  is  insulated  or  on  a 
line  of  symmetry,  the  zero  heat  flux  condition  holds  and  0  =  0. 


Convective  boundary 

Again  consider  the  boundary  illustrated  in  Figure  2.  When  surface  convection  is  present,  the 
heat  flow  at  the  boundary  is  Q  =  hAy(TA-Tt  j),  where  fA  is  the  ambient  temperature  outside  the 
grid,  and  h  is  the  surface  heat  transfer  coefficient.  The  heat  flow  between  node  /,/  and  the  sur¬ 
rounding  three  nodes  will  be  the  same  as  that  given  for  the  constant  flux  boundary.  The  steady- 
state  formulation  for  the  right-hand  boundary  with  Ax  =  Ay  -  As  is 


_1_  +  J_ 

*i-i.  i  *i,j 


(^i-t.i  -T.J  +  /— 


*i.M  *U 


(fu-l-T-j.j) 


1  I 

*1+1, j  *i,j 


(ri+t,j  -  ri,  f)  +  hAs(rA-Tu)  =  0. 


and  rearranging  gives, 


*i-t,i  *1.  j 


^JT  r"> ZlTX P'" 


*l.i-l  *i,i 


+  (_L_'+  J_  + 

\*l+1.i  *i,j/  \**-l,i  *»,i  *i,  j-I  *i,j 


+  hAs\  Tjj  =  -hAsT^ 


Again,  the  indices  may  be  suitably  rearranged  for  nodes  on  other  boundaries. 

Constant  temperature  boundary 

For  a  constant  temperature  node  on  any  boundary  or  inside  the  gride  we  have  T{  j  =  C,  where 
C  is  the  temperature  of  the  node. 

Semi-infinite  boundary 

The  semi-infinite  boundary  condition  represents  a  continuous,  uniform  material  extending  ad 
infinitum  in  one  direction,  with  a  known  temperature  at  infinity.  This  boundary  is  approximated 
in  the  finite  difference  computer  program  by  specifying  a  large  distance  between  the  boundary 
node  and  the  node  adjacent  to  it,  and  the  condition  requires  a  special  heat  balance  for  nodes  ad¬ 
jacent  to  the  boundary  as  well  as  for  the  boundary  nodes.  For  example,  consider  the  ri#it-hand 
semi-infinite  boundary  illustrated  in  Figure  3,  on  the  edge  of  a  grid  with  uniform  internodal 
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Figure  3.  Semi-infinite  boundary  approximation. 


distances.  The  temperature  at  infinity  is  represented  by  an  imaginary  node  /,  located  outside  the 
grid.  Dj-Ax  is  the  distance  from  the  internal  node  to  the  node  on  the  boundary,  and  also  the  dis¬ 
tance  from  the  boundary  node  to  the  imaginary  node  at  infinity.  As  previously  stated,  the  accur¬ 
acy  of  the  calculated  temperature  distribution  is  dependent  on  the  distance  between  nodes;  thus 
when  choosing  Dt,  we  should  choose  the  smallest  distance  which  may  still  be  considered  large  in 
terms  of  the  grid  size. 

The  heat  flow  between  node  i,j  and  each  of  the  surrounding  nodes  is  as  follows: 


(2Di-i)Ax-^r  (rM  r  rM) 


<?2 


Ay 


£>j  Ax 


f^i,  j-1  "  ^i,  j)> 


q3- 


(20,-1)  Ax  -  (7’j+i 


i-V’ 


Qa 


Ly'  D~Ax 


For  steady  state,  the  heat  balance  yields  (?,  +  Q2  +  Q3  +  Qa  =  0-  Allowing  Ax  =  A y,  we  find 
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Now  consider  the  heat  flow  for  the  square  node  adjacent  to  the  irregular  semi-infinite  node,  as 
illustrated  in  Figure  4. 
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Figure  4.  Node  adjacent  to  semi-infinite  boundary  node. 


The  heat  flow  between  the  node  and  each  of  its  adjacent  nodes  is 
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The  resulting  equation  for  steady  state,  for  Ax  =  Ay,  >s 
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(15) 


Constant  flux  corner 

Consider  the  corner  shown  in  Figure  5,  subject  to  the  constant  flux  condition  on  two  sides, 


Q3  =  'A  {4>y  +  0S)  As. 

<pj  is  the  heat  flux  per  unit  area  crossing  the  top  boundary,  0S  is  that  crossing  the  side.  For  steady 
state,  0,  +  (?2  +  Q3  ~  0.  In  general,  for  Ax  =  Ay  =  As, 


Figure  5.  Node  on  corner  of  finite  difference  grid. 


Convective  corner 

For  the  corner  shown  in  Figure  5,  subject  to  convection  on  both  sides,  convection  along  the 
top  surface  yields 

<?4=  h  f  (rA-rifj), 

and  along  the  right  side, 

<?3=  h  ^  (rA-ru). 
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With  (?,  and  Q2  the  same  as  for  the  previous  corner,  and  for  Ax  =  Ay  =  As,  the  equation  for 
steady  state  is  as  follows: 


Semi-infinite  corner 

Consider  the  node  pictured  in  Figure  6,  on  the  corner  between  two  semi-infinite  boundaries. 

F,  is  the  temperature  a  large  distance  to  the  right  of  the  grid,  and  is  the  temperature  away  from 
the  top  of  the  grid. 
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For  steady  state,  the  heat  flows  sum  to  zero.  Summing  the  above,  then  multiplying  each  term  by 
1/(2Dj-l)  and  allowing  A x  -  Ay,  g:ves  us  the  heat  balance  for  this  corner  node 


In  the  computer  program,  it  is  assumed  that  ^  '  ^i,j  r 

n 


-  V,-;;  -  %V-; 
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When  this  semi-infinite  condition  is  used  for  a  corner,  a  special  heat  balance  for  the  interior 
node  labeled  7$  in  Figure  6  must  be  used.  This  balance  will  be  the  same  as  that  given  by  cq  1  S, 
except  that  <?,  will  be  replaced  by 


<?i  = 


Then  the  full  heat  balance  is  given  by 
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Corner  with  one  side  constant  flux,  the  other  side  convective 

It  is  possible  to  have  a  corner  subject  to  convection  on  one  side  and  constant  flux  on  the  other 
side.  For  the  corner  in  Figure  5,  subject  to  constant  flux  on  the  right  side  and  convection  on  the 
surface, 


lhi-1'TJ’ 

(^i+1  ,j  “  ^i,j)  > 


Q 3  =  Vi  4>As  , 

C4=  h  y  (7-a-7-u). 


Let  Ax  =  Ay  =  As, 


+  Vi  hAs\  r, 


-  Yi  hAsTf,  -  Zi  0As  . 


(20) 


Corner  with  one  side  constant  flux,  the  other  side  semi- infinite 

As  illustrated  in  Figure  7,  let  the  top  right-hand  corner  of  the  grid  have  a  semi-infinite  boundary 
on  the  right  side  and  a  constant  flux  boundary  on  the  top.  Allow  the  heal  flow  per  unit  area 
crossing  the  top  of  the  grid  to  be  0.  Then  the  heat  flows  for  node  i,  j  are  given  by 


Q,  =  0(2 Dj  -  I)  Ax  , 
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Figure  7.  Node  on  a  corner  of  semi-infinite  boundary. 


The  total  heat  balance  for  steady  state  with  Ax  =  Ay  =  As  is  given  by 
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2(2Dr1) 

1+1,1  *i,i 


=  -  0  (2Dr1)  As 
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(21] 


As  is  the  case  with  the  other  semi-infinite  boundaries,  the  node  adjacent  to  the  semi-infinite 
node  must  be  given  special  consideration.  Again,  j>  is  the  heat  flux  per  unit  area  crossing  the  top 
of  the  grid. 


Allow  Ax  =  Ay  =  As.  The  heat  balance  equation  for  the  node  follows: 


/  mure  S.  Xode  udiaiaU  to  sc ■nd-infirdte  corner  node. 


IS 


The  heat  balance  for  that  node  is  given  by 


(24) 


Computer  program  development 

In  order  to  solve  a  particular  problem,  a  grid  whose  boundaries  coincide  with  those  in  the  prob¬ 
lem  must  be  set  up.  The  number  of  grid  points  used  is  largely  a  matter  of  trial  and  error;  as  pre¬ 
viously  stated,  the  smaller  the  As  between  the  nodes  (and  hence  the  more  dense  the  grid),  the 
more  accurate  the  solution.  Once  the  grid  is  set  up,  each  node  in  it  is  assigned  the  finite  difference 
equation  which  represents  the  heat  balance  between  that  node  and  its  adjacent  nodes.  The  resulting 
system  of  equations  may  then  be  solved  simultaneously  to  arrive  at  the  steady-state  temperature 
distribution  within  the  grid. 

SSCONDUCT,  the  program  developed  to  model  two-dimensional  steady-state  heat  conduction, 
is  made  up  of  four  parts:  1 )  a  data  gathering  subroutine,  2)  the  main  program,  3)  a  subroutine  to 
solve  the  system  of  equations,  and  4)  a  subroutine  to  locate  the  isotherms  in  the  resulting  tempera¬ 
ture  distribution.  The  grid  for  the  problem  is  represented  by  a  three-dimensional  array,  RAY  (I,  J,  K). 
The  first  two  dimensions  designate  the  spatial  location  in  Cartesian  coordinates.  The  third  dimen¬ 
sion,  K,  has  three  values.  RAY  (I, ),  1)  contains  the  temperatures  for  each  nodal  location  I,  j.  RAY 
(I,  J,  2)  describes  the  location  type  of  each  node.  Examples  of  location  types  include  a  node  on  a 
zero  flux  boundary,  a  constant  temperature  node,  a  variable  interior  node,  etc.  RAY  (I,  J,  3)  is  an 
index  of  the  material  type  of  each  node.  This  index  is  used  to  assign  conductivities  and  convection 
coefficients.  It  is  possible  to  have  a  different  material  at  each  point  in  the  grid. 


Data  subroutine 

Subroutine  SSDATA  was  written  to  provide  a  quick  and  easy  way  for  the  user  to  set  up  the 
grid,  without  having  to  worry  about  formats  in  the  data  file.  The  user  has  only  to  edit  SSDATA, 
following  directions  in  the  comment  statements  in  the  computer  program,  to  insert  the  desired 
values  of  the  variables.  When  SSDATA  is  run,  the  new  values  of  the  variables  are  put  into  the 
formatted  data  file  STSDAT.  The  main  program  then  reads  the  data  from  STSDAT. 

Main  program 

In  the  main  part  of  SSCONDUCT,  each  node  of  the  grid  is  examined.  The  conductivities  are 
figured  between  the  node  being  examined  and  the  adjacent  nodes,  and  the  coefficients  for  the 
node's  difference  equation  are  figured  and  stored.  A  grid  whose  dimensions  are  X  by  Y  contains 
XY  nodes;  hence  there  will  by  XY  equations  with  XY  unknowns.  Since  the  matrix  of  coefficients 
for  this  system  of  equations  is  banded,  with  all  entries  being  zero  outside  the  band,  only  the  entries 
of  the  band  need  to  be  stored.  The  bandwidth  is  2Y+1.  Thus,  instead  of  storing  an  XY  by  XY 
array  for  the  matrix  of  coefficients,  an  XY  by  2X+1  array  is  stored.  After  each  node  in  the  grid 
has  been  examined,  subroutine  BANDMX  is  called  to  solve  the  system  of  equations. 

Subroutine  BANDMX 

The  International  Mathematical  and  Statistical  Library’s  subroutine  LEQT1B  has  been  adopted 
for  solving  the  system  of  equations  which  has  been  stored  in  band  form.  Please  consult  Martin  and 
Williamson  (1967)  for  a  detailed  discussion  of  the  method. 
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Subroutine  tSOTHM 

This  subroutine  examines  the  final  temperature  distribution  in  the  grid,  performing  a  linear  in¬ 
terpolation  between  grid  point  temperatures  to  determine  the  spatial  coordinates  of  user-specified 
isotherms.  The  coordinates  are  written  into  file  POINTS  which  may  then  be  plotted. 


FINITE  ELEMENT  METHODS 
Introduction 

Finite  element  methods  arc  a  relatively  recent  addition  to  numerical  methods.  They  can  be  used 
to  obtain  approximate  solutions  to  governing  differential  equations  encountered  in  many  disciplines. 
They  were  first  used  by  the  aerospace  industry  during  the  1950's  for  structural  analysis  of  complex 
systems.  By  the  1%0’s  people  found  that  the  method  could  be  applied  to  a  broad  class  of  problems. 

Finite  element  models  use  a  piecewise  approximation  to  the  governing  differential  equations 
over  the  region  of  interest;  finite  difference  models  use  a  pointwise  approximation.  A  major  ad¬ 
vantage  of  the  finite  element  approximation  over  a  finite  difference  approximation  is  its  ability  to 
model  irregularly  shaped  boundaries  more  accurately.  The  si/c  of  the  elements  can  also  be  easily 
varied.  The  basic  shape  may  also  be  varied.  Figure  9  shows  how  an  irregularly  shaped  boundary 
would  be  modeled  with  both  the  finite  element  and  finite  difference  methods.  Each  of  the  models 
approximates  the  circular  boundary  ;  however,  the  finite  difference  model  uses  a  rather  crude  approx¬ 
imation  in  comparison.  Another  major  advantage  of  the  finite  clement  method  is  apparent  from 
Figure  9,  that  is,  we  can  easily  vary  the  element  size.  The  geometry  here  is  a  hall  symmetry  ol  a 
buried  pipe.  Wc  expect  the  temperature  gradients  to  be  much  greater  in  areas  nearer  the  pipe, 
and  by  varying  element  si/e  these  areas  can  be  modeled  to  any  degree  ot  accuracy ,  while  surround¬ 
ing  areas  with  smaller  temperature  gradients  need  not  be. 

The  finite  element  method  does  have  several  disadvantages.  Since  the  grid  is  usually  oddly 
shaped,  defining  all  the  coordinates  of  each  node  can  be  a  tedious  job.  Automatic -mesh-generation 


a.  Finite  difference  model. 


b.  Finite  element  model. 


Figure  9.  Finite  difference  and  Unite  element  grids  for  a  circular  boundary. 
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computer  programs  are  available  but  not  necessarily  easy  to  use  or  inexpensive.  Another  disadvan¬ 
tage  of  the  finite  element  method  is  its  complexity.  The  necessary  equations  cannot  always  be  en¬ 
tirely  based  on  physical  considerations,  as  can  finite  difference  equations.  Finite  element  computer 
programs  tend  to  be  complex  and  very  hard  to  modify.  In  spite  of  these  disadvantages,  the  finite 
element  method  is  widely  used  in  many  disciplines  of  engineering  today. 

Before  going  into  the  development  of  a  finite  element  model  for  heat  conduction,  we  will  first 
consider  the  basic  steps  which  must  be  followed  when  using  the  finite  element  method  to  solve 
any  problem. 

The  first  step  is  to  discretize  the  continuum,  th?'.  is,  divide  the  region  to  be  modeled  into  ele¬ 
ments.  Many  different  shapes  can  be  used  for  the  element  Here  we  will  consider  only  the  most 
simple  two-dimensional  element,  the  triangle.  The  nodes  and  elements  must  now  be  numbered. 

As  we  will  see  later,  this  must  be  done  carefully  in  order  to  get  the  most  economical  solution. 

The  next  step  is  to  select  the  type  of  interpolating  function  which  will  be  used  over  the  element 
to  represent  the  field  variable.  Normally,  polynomials  are  chosen  since  they  are  easily  differentiated 
and  integrated.  Since  we  will  be  using  a  simple  triangle  with  only  three  points  our  interpolation 
functions  will  be  linear. 

Next  we  need  to  find  the  equations  which  express  the  properties  of  each  element.  These  will 
be  in  matrix  form.  Several  approaches  are  possible.  We  will  rely  on  the  variational  principal,  which 
is  known  for  the  heat  conduction  problem.  This  approach  is  most  convenient,  but  in  some  cases 
a  variational  principal  does  not  exist  for  the  problem.  In  these  cases,  either  the  direct  approach, 
the  energy  balance  method,  or  die  method  of  weighted  residuals  (Galerkin's  appraoch)  may  be 
used  (Huebner  1975). 

Now  the  element  equations  can  be  assembled  into  the  global  equations  describing  the  entire 
region.  This  is  a  straightforward  procedure  easily  handled  by  the  computer.  The  system  of  equa¬ 
tions  must  be  modified  to  account  for  any  boundary  conditions  present.  The  result  is  a  system 
of  simultaneous  equations. 

The  next  step  is  to  solve  this  system  of  equations.  This  is  usually  the  most  time  consuming 
step  (computer  time).  With  linear  equations,  as  we  will  have,  many  methods  exist  for  solution. 

The  solution  is  somewhat  simplified  because  the  resulting  matrix  is  banded  and  svmetrical.  It  is 
also  positive  definite.  Round-off  errors  are  also  reduced  because  of  these  properties. 

Once  we  have  die  solution  we  may  want  to  compute  additional  quantities,  such  as  isotherms 
or  temperature  gradients  in  the  case  of  heat  conduction. 

In  the  next  section  we  will  discuss  the  procedure  used  at  each  of  these  steps  by  developing  a 
two-dimensional  finite  element  model  for  heat  conduction.  Various  types  of  boundary  conditions 
will  be  included  in  this  model. 

Application  to  two-dimensional  steady-state  heat  transfer 

In  this  section  we  will  develop  a  finite  element  model  for  heat  conduction  in  two  dimensions. 
The  simplest  two-dimensional  element,  the  triangle,  will  be  used  exclusively.  Linear  interpolation 
functions  will  be  used  within  the  elements  and  provisions  for  internal  heat  generation  will  be  in¬ 
cluded.  In  an  attempt  to  make  the  method  as  clear  as  possible  we  will  follow  the  sequence  of 
steps  outlined  in  the  last  section. 

The  first  step  is  to  divide  the  region  to  be  modeled  into  triangular  elements.  During  this  step 
we  should  try  to  concentrate  the  smaller  elements  in  the  areas  of  the  highest  gradients.  Larger 
elements  can  be  used  in  areas  of  smaller  gradients.  The  triangles  should  be  as  close  to  equilateral 
as  possible  in  order  to  promote  accuracy  of  the  solution.  As  an  example,  consider  a  problem  with 
a  geometry  similar  to  that  of  Figure  9.  The  region  can  be  divided  into  triangular  elements  in  a 
similar  fashion  as  shown  in  Figure  10. 

Now  the  nodes  need  to  be  numbered.  Care  must  be  taken  here  in  order  to  reduce  computational 
requirements.  As  stated  before,  the  simultaneous  equations  which  result  in  a  Finite  element 
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equation  .uc  banded  when  represented  in  matrix  lorm.  The  smaller  the  bandwidth,  the  easier  the 
equations  ate  to  solve  Storage  requirements  are  also  reduced  since  only  the  band  need  be  stored. 
The  numbering  ol  the  nodes  diiectly  allects  the  bandwidth  ol  the  matrix.  The  bandwidth  is  equal 
to  one  plus  the  largest  difference  between  node  numbers  lor  any  single  element  (out  of  all  the  ele¬ 
ments)  tor  a  problem  with  one  degree  of  freedom  (temperature)  at  each  node.  In  Figure  10  the 
nodes  have  been  numbered  in  a  way  that  we  hope  will  minimi/e  the  bandwidth. 

1  he  elements  must  also  be  numbered.  I  hey  may  be  numbered  in  an\  sequence,  however,  with¬ 
out  aflecling  the  computational  efficiency  of  the  solution.  The  computational  procedure  will  re¬ 
quire  that  we  know  which  node  numbers  are  associated  with  each  element.  We  must  also  be  able 
to  identity  what  material  a  particular  elements  is  in,  it  more  than  one  matetial  exists  in  the  problem 
A  table  '  t  the  format  ol  fable  I  will  accomplish  this  purpose. 


Table  1 .  Typical  element  data. 


Element 

Nodes 

Material 

number 

1 

2 

3 

type 

1 

to 

2 

1 

1 

2 

10 

11 

2 

1 

3 

11 

3 

2 

1 

4 

11 

12 

3 

1 

5 

12 

4 

3 

1 

• 

• 

• 

• 

• 

The  local  node  numbers  (1,  2,  3  in  Table  1 )  for  each  element  should  always  be  consistently  num¬ 
bered  with  a  particular  convention.  In  this  case  they  are  numbered  counterclockwise  around  the 
element 

We  also  need  to  define  the  spatial  coordinates  of  each  of  the  nodal  points.  A  simple  table  like 
Table  2  will  accomplish  this. 


Table  2.  Typical  nodal  point  data. 


Node 

X 

Y 

1 

0.0000 

-0.3280 

2 

0.1 256 

-0.3031 

3 

0.2319 

-0.2319 

4 

0.3031 

-0.1286 

5 

0.3280 

0.0000 

• 

• 

• 

The  coordinates  given  are  for  the  nodes  in  Figure  10.  The  origin  in  this  case  is  located  at  the  center 
of  the  pipe. 

Now  we  need  to  define  the  interpolation  functions  for  the  elements.  Certain  continuity  require¬ 
ments  must  be  met  by  the  element  type  and  the  interpolation  functions,  depending  on  the  type  of 
problem  being  studied.  It  is  beyond  the  scope  of  this  report  to  develop  these  requirements.  In  this 
case,  the  requirements  are  met  by  triangular  elements  and  the  linear  interpolation  functions  which 
will  be  used. 

The  interpolation  function  must  define  the  field  variable,  in  this  case  temperature,  over  the  en¬ 
tire  element.  That  is,  given  the  position  of  any  point  within  the  element,  the  element  interpolation 
function  gives  an  estimate  of  the  temperature  at  that  point  The  form  of  the  linear  interpolating 
polynomial  is 

T=a1  +  a2x  +  a  3y  (25) 

where  T  is  the  temperature,  x  and y  are  the  coordinates  and  a,,  a2  and  are  constants. 

At  local  nodes  1,  2  and  3  of  die  element  the  following  conditions  must  be  satisfied  by  the  inter¬ 
polating  polynomial: 


/  /  (  at  v  A  ,  and  y  -  V' 3  .  (26) 

Using  each  ol  these  conditions  and  solving  lor  /  at  each  node  gives  us  the  following  set  ot  equations: 


-  a 

+  n.A, 

/ , 

a 

*  ^2  ^2  +  V  2 

/ 

'  } 

a 

i  ^  (J  i  A  ^  +  qa  V  ^  . 

(2?) 

Now  we  can  solve  this  set  ol  equations  tor  the  constants  o  | ,  a,,  a  3,  yielding 

al 

~  |(V_,V3  -  \  ,  >  ,)/ ,  MV(V,  -.V,>  ,) /,  +  (A  ,  V  -  A  ,  >  ,)/,! 

a. 

.  J 
2 

y  l(A  ,>■,)/,  MA(>  , )/  s  *  (V,  ->_,)/,! 

~  |(\,.v,)/,  +  (A  ,  A  t)  /  j  +  (A  j  -  V,)/.|  . 

(2b) 

where  1 

is  the  area  ol  the  triangular  element  and  2  1  is  given  In  the  deleiminant  ol 

the  following 

3s  3  matrix: 

1  A,  V  ,  j 

2-1 

1  Vs  >  s  | 

'  Vt  >  s  j 

(2‘» 

Now  we  need  to  find  the  inlet polation  (unction  in  terms  ol  the  temperatures  at  each  node  and 
the  element  topographs .  We  can  do  this  by  substituting  out  expression  lor  the  constants  q  j  ,  o_,, 
a,  into  the  otigin.il  interpolation  function  (eq  25).  Ihe  result  will  be  an  equation  lot  / 


/  If  I 

V,  /  ,  *  +  v,/t 

Ol 

/<•> 

ivi  M:pl 

(30) 

Hete  the 

superscript  (e)  means  "within  the  element"  and  ,V| 

,  V,  and  A"3  ate  the  shape  (unctions; 

thev  ate  given  bv  the  lollovving  equations: 

v, 

l,; t  *  ^t '  +  ‘,y| 

V, 

K>  •  •>_>'  *  *  d  l 

v.t 

I"',  +  /?t'  f  ‘  c'  l  • 

(31) 

1  he  constants  in  these  equations  ate  given  by 

"l 

V:>  r  -  \  ,V, 

a  s 

V  s'  ,  -  A , V  , 

"t 

A  ^  2  '  A  2  ^  \ 
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C ,  -  x3  -  x2 

c2=  A,  -  X3 

ci=X2-X].  (32) 


It  is  easily  shown  by  substituting  bad.  into  cq  3)  that,  lor  any  particular  node,  the  corresponding 
shape  function  takes  on  a  value  of  unity  at  that  node  and  zero  at  the  other  two  nodes  and  the  line 
connecting  them. 

Now  we  can  use  the  information  in  Tables  I  and  2  to  assemble  the  individual  element  equations 
into  a  global  system  of  equations.  The  result  is  a  set  ol  piecewise  continuous  equations  to  approxi¬ 
mate  temperature  over  the  entire  region.  Our  next  step  is  to  relate  this  to  the  physical  problem  ol 
heal  conduction  in  a  solid  medium. 

Several  methods  are  available  for  finding  the  element  equations,  as  mentioned  earlier.  I  ot  this 
problem  we  will  use  the  variational  principal,  which  is  known.  A  discussion  ol  the  calculus  of 
variations  is  beyond  the  scope  of  this  report.  The  interested  reader  is  referred  to  Pais  ( l%2), 
Schechter  (1967),  Huebncr  (1975),  Segerlind  (1976)  and  Zienkiewic/  (1977). 

The  calculus  of  variations  simply  states  that,  if  wc  can  lind  the  so-called  '‘functional"  which 
corresponds  to  the  governing  differential  equation  and  boundary  conditions  ol  the  problem,  the 
minimization  of  that  functional  requires  that  its  corresponding  diflcrential  equation  and  bound¬ 
ary  conditions  be  satisfied.  For  isotropic  two-dimensional  steady-state  heal  transfer,  the  governing 
differential  equation  is 


(-») 


where  Q(  is  the  internal  heat  generation  within  the  region.  Constant  temperature  boundary  con¬ 
ditions  are  simply 


/  7  on  5 1 

where 

Ty  =  constant  temperature  of  the  boundary  5, 

5,  =  segment  of  boundary  at  lg. 

Constant  heat  flux  boundary  conditions  arc  given  by 
A  !^*+A’  ^V'^0  »"  *2 

where 

lfx  and  t!v  =  direction  cosines  of  a  vector  perpendicular  to  .52. 
S2  -  surface  subject  to  the  constant  heat  flux  0. 

Finally,  the  conditions  on  a  convective  boundary  arc 

A  dT^+A  ^  53 

where  S3  is  the  boundary  segment  subject  to  convection. 


(35) 


(36) 
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The  entire  boundary  of  the  region  is  made  up  completely  ot  segments  S,,S2  and  S3.  Radiative 
boundary  conditions  are  not  considered  here. 

Now  we  need  to  use  the  variational  principal  to  find  the  element  equations.  The  functional 
matching  eq  33  with  boundary  conditions  (eq  34,  35  and  36)  is  (Schechter  1967,  Huebner  1975, 
Segerlind  1976) 

'(T)  =  \  /  [*(|92  +*(t£)2  -2Qj\dxdy  *  I  |0T+l/2Mr-rA)2US4  (37) 

A  S4 

where  54  is  the  union  of  surfaces  52  andS3. 

We  can  define  T over  the  element  by  recalling  eq  30 

7ie>  =  \N\  jr{te>  . 

To  minimize  the  functional  /(T),  which  will  in  turn  satisfy  governing  differential  equation  and 
boundary  conditions,  we  can  set  its  derivative  with  respect  to  temperature  /  equal  to  zero.  Be¬ 
cause  our  interpolation  functions  meet  the  continuity  requirements  as  discussed  earlier,  the  total 
integral  /( T)  is  the  sum  of  the  integral  l(T^c^)  for  all  the  elements.  If  we  have  M  elements, 

M 

t(T)  =  £  /(7l<'>).  (38) 

e  =  I 

Because  the  temperatures  at  each  node  of  an  element  are  independent,  the  partial  derivative  ol 
the  integral  l(T (*'*)  with  respect  to  each  nodal  temperature  must  be  zero 


a /  (r*‘->)  _  a/(r<<'>)  _  a/ (/<*•> 
ar,  ‘  a/2  ar3 


If  node  I  is  on  the  boundary  of  surface  S4,  we  have 


a/(  rU'>)  _ 


ar, 


arfH 

a 

/ar<^\ 

,  a /■<*•) 

+  t  — 

a 

t9— i 

*  3a 

ar, 

Ids  / 

3> 

a/', 

*  a.  / 

^  ar(f *  |  f  /  ar<(,) 

+  <?,  — — I  <^(<)  +  /  (o  tt- 

s>>  ’ 


a/, 


(40) 


Similar  equations  apply  to  other  nodes  on  boundary  54.  for  nodes  not  on  boundary  S4,  the 
integral  over  54  is  dropped.  The  derivatives  in  the  above  equation  can  be  evaluated  by  considering 
eq  30.  After  substituting  these  into  eq  40,  wc  find  the  result  tor  node  I  on  boundary  S4  is 


a  /  (rH) 
ar, 


\r\ 


a/v, 

by 


+  Q.N, 


dA  h'l 


+ 


/ 


*4 


+/»|tv|  /v,  jr-rA}|  ds4(f»  . 


(41) 
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The  element  equation  for  each  element  is  now  found  by  combining  the  resulting  three  nodal 
equations, 


a /  |(‘>  _  \'d)i(T(e>)  a/(r(f>)  a^r^i)' 

an*’) 


L  ar,  •  a r2  *  a r3 


!] 


transpose 


=  [/c|(4')  +  {/?}(f>  +  (*s4](f)  M,f) 

=  (*,<0  +  [Ks4]<<’j  {rf<f>  *  }«}(*>  =o. 


(42 


Both  the  (AC  ( ( e )  and  |ACj4  |  (e)  are  3x  3  matrices  whilejAfjl‘4  and  {/{(‘4  arc  3*  1  vectors.  Typical 
terms  within  them  arc 


c  /3  N.  a  Nj  a  N.  3AS\ 

/  *(aT  at2  ♦  TT  x) 

A<‘l 

R>  =  f  (?,A2  (//>(<’)  +  /"  0.\\c/54<''> 

A<f>  ‘  (,<'> 

fs.  =  f  /»A,Al£y54(*l  . 

12  "c.(*’f 


(43) 


Now  each  of  these  equations  tor  an  element  can  he  assembled  into  the  global  system  ol  equa¬ 
tions.  The  procedure  for  doing  so  is  straightforward.  If  we  have  n  nodes  in  out  problem,  out  re- 
suiting  matrix  will  have  dimensions  of  n\ n.  Using  the  data  in  Table  1  we  simply  add  all  ol  the 
enti  ics  ol  each  of  the  element  matrices  info  the  corresponding  global  position  in  the  nxo  sy  stem 
rriJir ix.  The  same  basic  procedure  applies  to  finding  the  global  R  vector. 

The  assembled  system  of  equations  needs  to  be  modified  to  account  for  constant  temperatuie 
boundary  conditions.  Although  exact  methods  are  available  for  this  proceduie,  we  will  use  an 
approximate  technique  which  is  much  simpler  to  program  for  computer  execution.  Befoie  modi¬ 
fication,  eq  42  is  rewritten  into  the  form  tor  solution 

IM  {/MM  (44) 

where  {/,{  is  now  the  sum  of  the  global  jA>  (  and  |Ay  ,  |  matrices  and  {/  |  is  simply  -  |A.’|  It  sse 
have  a  constant  temperature  boundary  condition  at  a  given  node  we  simply  multiply  the  corre¬ 
sponding  diagonal  term  in  the  |A  j  matrix  by  a  relatively  large  numbet  and  also  replace  the  corre¬ 
sponding  term  in  the  {/'{vector  by  the  product  of  the  boundary  temperature  and  this  modified 
diagonal  term.  In  this  case  we  have  used  1 0 1  ^  as  fhe  large  number.  I  his  is  an  approximate  pio- 
ceduie  but  it  will  yield  good  results  as  long  as  the  boundary  temperatuie  specified  is  not  very  small 
(Segerlind  1976). 

Now  we  are  ready  to  solve  for  the  unknown  temperatuies.  II  our  problem  has n  nodes  our  |A,  | 
matrix  will  be  an  nxn  matrix  and  both  j/j<ind{/  {will  bens  I  vectors.  It  becomes  obvious  that 
for  large  problems,  often  having  hundreds  of  nodes,  this  is  the  most  time  consuming  step,  f  or 
this  reason  it  is  imperative  that  our  solution  method  be  efficient.  The  fact  that  our  resulting  |A  ] 
matrix  has  some  special  characteristics,  as  mentioned  earlier,  greatly  simplifies  and  speeds  the 
solution. 

Many  techniques  are  available  tor  solving  linear  setsol  equations  such  as  we  have  here.  I  he 
interested  reader  may  refer  to  one  of  the  many  texts  on  the  subject,  including  I  orsythc  et  al. 
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(1977)  and  Gerald  (1978).  The  method  used  here  is  in  the  form  of  a  computer  program  subroutine. 
The  subroutine  is  called  MCHB  and  is  part  of  the  IBM  SSP  (Scientific  Subroutine  Package),  which 
is  available  on  many  computer  systems.  The  matrix  which  is  provided  to  MCHB  must  be  in  com¬ 
pressed  form,  consisting  of  only  the  main  diagonal  and  the  upper  portion  of  the  band,  which  is 
symmetrical.  The  terms  are  stored  in  rows  in  successive  storage  locations.  It  returns  the  tempera¬ 
ture,  after  solution,  in  the  same  space  that  the  jcf  vector  was  supplied  in.  The  Cholcsky  decompo¬ 
sition  technique  is  used  by  MCHB  to  solve  the  matrix  equations. 

All  the  steps  of  solution  for  a  two-dimensional  steady -state  heat  transfer  problem  using  the 
finite  element  method  have  now  been  outlined.  The  only  remaining  step  is  to  compute  additional 
quantities  based  on  the  resulting  temperature  distribution.  Olten,  the  location  of  a  particular  iso¬ 
therm  rru\  be  of  interest.  A  simply  method  of  finding  the  approximate  location  of  an  isotherm  is 
simply  to  interpolate  linearly  along  the  element  boundary  between  nodes  whose  temperatures 
bound  the  isotherm  of  interest. 

In  the  next  section  we  will  discuss  the  computer  program  based  on  the  theory  developed  in 
this  section. 

Finite  element  computer  program 

A  finite  element  computer  program,  called  FEHEAT,  has  been  developed  to  model  two-dimen¬ 
sional  steady-stale  heat  conduction.  This  program  is  based  in  part  on  a  program  presented  by 
Huebner  ( 1975).  The  following  types  of  boundary  conditions  have  been  included: 

1.  Constant  temperature 

2.  Specified  heat  flux  (including,  ot  course,  zero  flux  insulated  boundaries) 

3.  Convective  heat  transfer. 

I  EHEAT  also  has  the  capability  to  account  for  internal  heat  generation. 

The  program  consists  ol  a  main  program  and  live  subroutines.  The  main  program  is  responsible 
for  the  following  operations'. 

1.  Providing  dimensioned  and  common  storage  spate  lor  the  necessary  matrices,  vectors,  and 
scalars 

2.  Initializing  all  of  the  above  to  zero,  as  required 

3.  Determining  the  size  of  the  problem  and  type  of  boundary  conditions  present 

4.  Reading  in  all  input  data  and  printing  them  out  tor  confirmation. 

5.  Determining  the  bandwidth  of  the  |Al  matrix 

6.  Calling  the  appropriate  subroutine  for  solution  and  isotherm  location 

7.  Printing  out  resulting  temperatures’  and  isotherms'  coordinates 
The  names  and  purposes  of  the  five  subroutines  are: 

1.  ISM  This  subroutine  forms  the  \K  J  b'l  matrix  (element  “thermal  stiffness”  matrix).  A 
different  type  ol  material  may  be  used  in  each  clement.  TSM  computes  the  constants  in  the  linear 
interpolation  lunciions  as  well  as  the  element  area.  The  “thermal  stiff  ness  matrix  is  modified 

to  account  lor  convection  and  the  internal  heat  generation  at  each  node  is  computed  here  for  use 
in  subroutine  FRHS. 

2.  FRHS  This  subroutine  forms  the  jF|vcctor  (right  hand  side).  It  begins  with  the  contribu¬ 
tion  of  internal  heat  generation  at  the  affected  nodes  and  adds  to  that  the  effect  of  specified  heat 
fluxes  and  convective  boundary  segments.  It  is  then  modified  for  the  constant  temperature  bound¬ 
ary  conditions  by  the  procedures  outlined  in  the  previous  section. 

3.  FORMK  This  subroutine  assembles  the  clement  matrices  to  form  the  global  |Aj 

matrix.  Constant  temperature  boundary  conditions  arc  accounted  for  by  modification  of  the 
corresponding  diagonal  element  as  outlined  earlier.  The  resulting  matrix  is  stored,  in  rows  first 
then  in  columns,  in  consecutive  storage  locations,  with  only  the  diagonal  and  upper  portion  of 
the  band  being  stored,  as  required  by  MCHB.  This  procedure  greatly  reduces  storage  requirements 
for  banded  matrices. 


4.  MCHB-This  subroutine  solves  the  system  of  equations.  The  (A'|  matrix  must  be  stored  in  the 
form  explained  earlier.  Double  precision  is  used  on  several  crucial  quantities.  This  subroutine  is 
part  of  the  IBM  Scientific  Subroutine  Package  available  on  many  computers. 

5.  ISOTHM-This  subroutine  finds  the  location  of  any  number  of  specified  isotherms  at  any 
temperature.  It  examines  each  clement  and  performs  linear  interpolation  between  adjacent  nodes 
within  an  element  which  bounds  the  temperature  of  interest.  Its  results  arc  sets  of  coordinates  for 
each  temperature  specified  by  the  user.  These  results  can  be  plotted  with  a  simple  plotting  program. 

More  will  be  said  on  the  computer  program  and  the  preparation  of  input  data  in  the  next  section 
where  we  examine  its  application  and  results  for  a  classical  problem.  A  listing  of  the  computer  pro¬ 
gram,  including  all  the  subroutines  and  sample  input  files  and  output  data,  is  given  in  Appendix  B. 


PROGRAM  VERIFICATIONS  AND  COMPARISONS 

Each  computer  program  was  run  modeling  two  steady-stale  two-dimensional  problems:  the 
case  of  a  rectangular  plate  subject  to  a  uniform  temperature  on  three  sides  and  a  sinusoidal  temper¬ 
ature  distribution  across  the  fourth  side,  and  the  case  ol  a  buried  pipe.  Analy  tical  solutions  exist 
for  both  of  these  problems.  In  the  Rectangular  plate  section,  the  rectangular  plate  problem  is 
solved  analytically,  and  the  results  are  compared  to  the  results  of  each  program.  In  the  Huried 
pipe  section,  the  analytical  solution  of  the  buried  pipe  is  presented,  along  with  comparisons  of  the 
computer -genetated  solutions.  In  the  Semi-infinite  condition  section,  a  one-dimensional  problem 
is  solved  by  the  finite  difference  program  to  illustrate  the  use  ol  the  semi-infinite  boundary. 

Rectangular  plate 

Consider  a  rectangular  plate  that  has  three  ol  its  sides  at  a  lived  lempetature  and  has  a  sinu¬ 
soidal  temperature  distribution  across  the  fourth  side.  The  following  well-known  analytical  solu¬ 
tion  exists  to  find  the  steady-state  temperature  distribution  within  the  plate. 

Assume  the  plate  is  sufficiently  thick  so  that  the  end  el  feels  may  be  neglected.  1  he  thermal 
conductivity  is  constant  throughout  the  plate.  We  examine  a  cross  section  of  the  plate,  and  use  a 
shitted  temperature  0  -  1  -  1 0.  The  heat  flow  obeys  Laplace’s  equation 

If®  ♦  O  =  0 

Tv-  cV-’ 

The  boundary  condition  on  the  lop  surface  is  0  -  0m  sin  rr\  /If;  II  -  0  lot  the  three  other  sut- 
faces.  This  gives  us  the  following  boundary  conditions: 

(0  0(O.y)  =  O 

(2)  0(M',y)  =  0 

(3)  0  (a  ,  0)  =  0 

(4)  0  (a,  //)  =  0m  sin  ~  .  I  »5) 

We  begin  by  assuming  a  solution  of  the  form 

0(x,  y) X(x)Y{y)  . 

Then,  from  Laplace’s  equation 

L  *  1  b2Y 

~  *  bx2  ~  y  by2 
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Wc  then  proceed  as  usual  with  the  method  of  separation  of  variables,  where  the  separation  constant 
is  X* : 


d-’A 


v  a* 

J-.X 
c/s  - 


X-, 


3 ~v 

dv~ 


\2 


+  A-.Y  =  0,  —  -  A*'K  =  0 

dy 2 


A'  -  C  |  cos  As  +  C->  sin  As,  V  -  C  ,e"A--'  +  C, e*1 

0  ■=  (C ,  errs  Aa  +C,  sin  As)  (C3e-X‘  +  C4cM  )  . 

Apply  the  boundary  conditions,  from  eq  4 5 1 , 

C,  (C  ,e-X'  +  C.,ex‘ )  -  0 
C,  -  0  . 
from  eq  45  >, 

(C\  cos  Av  +  C,  sin  As)  (C,  +  C, )  -  0 

-^4 

f  rom  eq  45  ,,  with  results  of  eq  45,  and  45,, 
sin  Alt  (eM  -e'M  )  ■'  0 
2C\ ,  sin  Alt' ->in It  Ay  =  0  . 

Ihis  requires  that  sin  Alt'  '  0  or  A  ~  tin) if  (/;  is  an  integer).  Because  ol  the  linearity  of  Laplace's 
equation,  0  can  he  written  as  a  sum  of  an  infinite  series 


"  £  f-„ 

n  » 


mr\  .  n  ttv 
sm  —  smh  -jjr- 


where  the  constants  have  been  combined. 

Now  apply  the  boundary  condition  from  eq  45.,: 


,  rrs  v  nn s  .  ,  rrrr// 

"m  Mn  jV  Cn  sin  —  smh 

n  :  I 
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This  hold-,  only  it 
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and  if  C2  =  C3  =  ...  =  Cn  =  0,  then  the  final  temperature  distribution  in  the  plate  is  given  by 


sinh  SL 

„  ..  _ W  .  nx 

0  0^  ~~Hh  sm  ¥ 

sinh  — 

W 


sinh 


T  =  r„  +  f „ 


sinh 


ny 

W  rrx 

Vh  sm  ¥ 


(46) 


For  the  computer  comparisons,  let  7 0  =  I00°C,  /m  =  50°C,  H  =  8  m,  and  M'  =  12  m.  The  results 
were  calculated  tor  every  1  m  increment  of  space  in  the  rectangle,  and  the  results  are  shown  in 
Table  3. 


Table  3.  The  rectangular  plate  solution. 

0*  7  2  .?  1  5  ft  7  X  9  II)  II  U 


0  100.00*  100.00  100.00  100.00  100,00  100.00 

1  100.00  100.8(,  101.66  102.34  102.87  103.20 

2  100.00  101.77  i  03.43  104.84  10S.93  106.62 

3  100.00  102.81  105.43  107.68  1 09.41  110.49 

4  100.00  104.04  107.81  111, OS  113.53  115.09 

5  100.00  105.55  110.7  3  1  15.17  118.58  120.73 

6  100.00  107.45  114.39  120.35  124.92  127.80 

7  100.00  109.85  119.04  1 26.92  132.97  136.78 

8  100.00  112.94  1  25.00  1  35.36  143.30  148.30 


100.00  100.00  100.00  100.00  100.00  100.60  100.00 

10.3. .31  103.20  102.87  102.34  101.66  100.86  100.00 

106.85  106.62  105.9.3  104.84  103.4  3  101.77  100.00 

110.86  110.49  109.41  107,68  105.43  102.81  100.00 

115.62  115.09  113.53  111.05  107.81  104.04  100.00 

121.46  120.7.3  118.58  115.17  I  1 0.73  105.55  100.00 

128.78  127.80  124.92  120.35  114,39  107.45  100.00 

138.08  1.36.78  132.97  12 6.92  119.04  109.85  100.00 

150.00  148.30  14.3.30  135.36  125,00  112.94  100.00 


•Increments  in  meteis. 
♦Values  in  'C. 


This  problem  was  run  on  SSCONDUCT,  the  finite  ditfcrence  program,  using  a  1 7*  25  grid  (425 
nodes)  with  an  inter  nodal  spacing  of  0.25  m.  The  results  compare  to  within  0.01  °C  ol  the  analyti¬ 
cal  solution. 

f  F.HEAT,  the  finite  element  program,  also  used  425  nodes.  The  elements  weie  isosceles  right 
triangles  with  sides  0.25  m  long.  The  total  number  of  elements  used  was  768.  The  results  from  this 
model  predicted  the  tempeiature  to  within  0.02  C  ul  the  analytical  solution  throughout  the  region. 

Buried  pipe 

The  solution  to  the  steady-state  problem  of  a  buried  pipe  may  be  lound  from  vector  lield  theory 
by  superposition  of  the  potentials  of  two  line  sources  with  equal  and  opposite  strength. 

In  general,  the  gradient  of  a  scalar  field  is  a  vector.  In  the  case  of  a  line  soutco,  we  assume  that 
the  source  is  of  uniform  strength  all  along  its  length.  Then  the  radial  vectors  arc  those  which  are 
of  interest.  In  the  case  of  heat  flow,  allow  the  radial  vectors  to  be  denoted  by  q.  The  scalar  poten¬ 
tial  field  is  the  temperature  distribution 

(47) 

where  k  is  the  thermal  conductivity,  and  q  is  given  by  q  -  WI2vr(r  ) ,  where  W  is  the  source  strength 
per  unit  length  (W/m),  r  is  the  distance  from  the  line  source  to  the  point  under  consideration,  and 
r  is  the  unit  vector  in  the  radial  direction  pointing  from  the  line  source  towards  the  point  under 
consideration,  as  shown  in  Figure  1 1. 
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Figure  12.  Superposition  ot  two  line  sources  (•  -heat  source  below  the  ground;  heat 
source  above  the  ground). 

Upon  substitution  into  eq  47: 


w  -  a  / 

2  nr  —  hr 


r 


mF  r«'- 


=  7-ro  =  0. 


(48) 


Now  consider  the  situation  illustrated  in  F  igure  12.  IF  is  the  heat  source  strength  below  the 
ground;  we  also  imagine  a  source  of  equal  and  opposite  strength  an  equal  distance  above  the 
ground.  The  field  lines  are  the  dotted  lines,  the  isotherms  are  represented  by  continuous  lines. 
For  the  bottom  source 
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and,  for  the  top  source, 

,  _  W  r2 

^  ~  ™  'n  77  ' 

The  resultant  field  (temperature  distribution)  may  be  obtained  by  superposing  the  potentials  Irom 
these  two  fields, 


W  i  r2  r\  \ 

0'  +0^=  2nk  (lnr  -|nr) 


.  W  .  r2 
<>  -  r —  In  — 
2nk  r , 


/  -  T. 


2  rtk 


(49) 


Isotherms  are  lines  of  constant  0-  When  /  _>  /  , ,  the  straight  line  between  the  two  I  ields  is  given  by 

W 

0  -  r-j-  In  I  =  0 
2nk 

r- td  --  0 
1  -  rn . 


Thus  the  reference  temperature  is  given  at  the  ground  surface. 

Now  put  the  problem  on  a  system  of  Cartesian  coordinates,  as  illustrated  in  I  igurc  I  f. 


t' 


/  iQiiw  13.  Cartesian  inordinate*  tot 
sour  1  e  problem. 


The  source  of  strength  U’  is  located  at  the  point  (t),  -(/)  I  torn  I  igure  14,  the  distances r2  and 
r,  arc  given  by 

r7  +  U'-V|)2 


}() 


Substitute  these  results  into  eq  49 


T~T0 


IT  v^i'  + 


4rr*(T-T0)  _  xf  +  (d-y,)2 

w  *2  +  (d+yt)2 


e\p 


/4  nk 
\  W 


*i  +[d-y,)2 

*i  (tz-t-xi  )- 


(50) 


This  gives  a  general  equation  for  the  temperature  at  any  point,  given  the  ground  surface  temper¬ 
ature  ro,  the  conductivity,  and  the  source  strength.  It  will  be  shown  below  that  the  resulting  iso¬ 
therms  have  a  circular  shape.  We  will  let  the  pipe  be  represented  by  one  of  the  isotherms.  Given 
the  temperature  of  that  isotherm,  we  may  calculate  the  resulting  temperature  distribution,  with  the 
temperature  of  the  pipe  and  the  temperature  ol  the  ground  surface  also  given. 

For  a  given  isotherm,  the  left-hand  side  of  cq  50  is  constant.  Let 


then 


r  = 


'i  +(r/-y,)' 

'l  +  (c )2 


After  the  fraction  is  cleared  in  this  equation,  we  complete  the  square  to  arrive  at 


■v?  * 


-1  • 


Thus  each  isotherm  is  a  circle  with  center  c  and  radius  r,  as  given  below. 

<■  =  («,  ^-rf) 


r=  2 d 


vL 

c-l 


For  comparison  with  the  computer  models,  consider  the  problem  of  a  pipe  with  a  0.1-m  radius  at 
100°C  buried  at  a  depth  of  1  m  (center  of  pipe);  the  ground  surface  is  at  10°C. 

The  pipe  is  the  1 00°C  isotherm  whose  center  is  at  1  =  d  (cQ  + 1  /c0  - 1 ) 


r  =  0. 1  =  2d 


co'] 
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Substitute  for  d  to  obtain 


0.1 


sfc~o 

(c0-1) 


cl  -  398  cD  +  1  =  0 


cQ  397.9975 
d  =  —— — —  =  0.995  . 

C„  +  1 


Then 

397.9975  =  exp|4n*/W'(100-10)  | 
5.986  =  360  nk/W 
W  =  360  rr/t/5.986  . 


For  any  isotherm, 

C  =  exp  (4rrA(7"-7"0/(T)  |  =  exp|4rr£(7  -  ro)  (5.986)/360  nfc ) 

C=  exp  |0.0665  (T-  ro)|  . 

Thus,  Riven  a  point  (a,  y)  of  interest,  its  temperature  may  be  calculated  from  the  lollowinR  equa¬ 
tion; 


t.0.0665(7'-l  0) 


x\  +  (0.995 -y,)2 
*?  +  (0.995  +  >'|  )2 


a?  +  (0.995  -y,  )2 

T=  10+ 15.033%  In  — — —  .  (51) 

a  l  +  (0.995  +y,)2 


The  region  modeled  by  the  computer  programs  was  a  rectangular  area  extending  from  the  giound 
surface  to  3  m  below  the  surface,  and  2  m  to  each  side  of  the  pipe.  Because  there  is  a  vertical  line 
of  symmetrv  •  Mending  from  the  ground  surface  down  through  the  center  of  the  pipe,  the  com¬ 
puter  progra:  modeled  only  half  of  the  area,  and  used  a  zero  flux  boundary  along  the  right  side, 
which  passes  through  the  center  of  the  pipe. 

The  finite  difference  program  used  a  31  x  21  grid  (651  nodes),  with  an  internodal  distance  of 
0. 1  meters.  The  pipe  was  represented  by  four  constant  temperature  nodes.  The  left  side  and 
bottom  boundaries  were  assigned  the  semi-infinite  boundary  condition.  See  Appendix  A  for  the 
input  and  output  for  the  program;  the  grid  is  labeled  and  printed  in  the  output  file,  STDOUT. 

This  problem  required  88  CPU-seconds  (Central  Processing  Unit)  on  CRREL’s  PRIME  400  com¬ 
puter,  a  cost  of  $3.08. 

The  finite  element  program  used  104  nodes  ( 166  elements).  The  elements  were  triangles  with 
sides  of  varying  lengths,  as  shown  in  Figure  14.  The  left  side  and  bottom  boundaries  were  assigned 
a  constant  temperature  boundary  condition.  The  finite  element  problem  required  26  CPU-seconds 
on  the  same  computer,  and  cost  $0.91. 

The  three  solutions-  analytical,  finite  difference,  and  finite  clement  are  graphed  in  Figure  1 5 
for  the  30°,  50°  and  70°C  isotherms.  Also  shown  arc  the  finite  difference  and  finite  element 
100°C  nodes,  which  represent  the  pipe.  Note  that  the  finite  clement  method  is  better  able  to 
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Figure  15.  Comparison  of  finite  difference  (*)  and  finite  ele¬ 
ment  (  )  results  for  buried  pipe  problem  (the  solid  line  is  the 
ana/vtical  solution). 


model  the  circular  surface  of  the  pipe,  while  the  finite  difference  method  used  only  a  four-node 
representation.  Both  methods  represent  the  temperatures  above  the  pipe  very  well,  but  the  finite 
difference  method  is  slightly  more  accurate.  For  regions  to  the  left  of  the  pipe  and  below  the  pipe, 
there  are  inaccuracies  in  both  of  the  computed  solutions,  mainly  due  to  the  effect  of  the  semi¬ 
infinite  boundaries.  The  finite  difference  solution  overpredicts  the  distance  of  the  various  isotherms 
from  the  pipe,  but  is  the  more  accurate  solution  in  this  case.  The  finite  clement  method  underpre¬ 
dicts  the  isotherm  locations.  The  internodal  spacing  used  could  be  reduced  in  both  cases  should  a 
more  accurate  solution  be  necessary. 

The  semi-infinite  condition 

The  semi-infinite  boundary  condition  in  the  finite  difference  program  is  modeled  by  use  of  a 
long  rectangular  node  at  the  boundary.  The  temperature  at  infinity  is  known.  To  verify  this 
approach,  consider  the  one-dimensional  problem  of  a  semi-infinite  plane.  The  temperature  at  the 
edge  is  !15°C,  and  the  temperature  at  infinity  is  0°C.  If  we  allow  "infinity"  to  be  a  large  but 
finite  distance,  then  a  linear  temperature  distribution  between  0°  and  1 1  5°C  would  be  expected 
at  steady  state. 

A  grid  17  nodes  long  was  used  to  model  this  situation  (Fig.  16).  There  are  16  nodes  of  rcgulai 
shape  in  the  grid,  and  one  semi-infinite  boundary  node.  Let  OS  be  the  internodal  distance  for  the 
interior  of  the  grid,  and  100-Z3S  be  the  distance  from  the  last  square  node  to  the  node  at  infinity. 
Listed  in  Table  4  are  the  expected  temperatures  for  each  node  and  those  calculated  by  SSCONDUCT. 


Tabic  4.  One-dimensional  semi¬ 
infinite  comparison. 
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figure  16.  Semi-infi¬ 
nite  verification  problem. 
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The  error  at  node  1 7  is  small,  considering  that  the  distance  modeled  by  the  one  rectangular  node 
is  almost  seven  times  the  size  of  the  rest  of  the  grid.  Errors  encountered  when  using  the  semi-infin¬ 
ite  condition  in  two  dimensions,  however,  may  be  much  larger,  as  may  be  seen  from  the  buried 
pipe  example.  In  general  the  semi-infinite  condition  should  be  used  in  two-dimensional  problems 
only  in  regions  where  the  temperature  gradient  is  small  and  precise  knowledge  of  the  temperature 
distribution  is  not  critical. 


INSTRUCTIONS  FOR  USE  OF  COMPUTER  PROGRAMS 

This  section  is  provided  for  those  who  wish  to  use  either  one  ol  the  two  heat  conduction  pro¬ 
grams.  SSCONDUCT  is  the  finite  difference  program,  and  FEHf.AT  is  the  finite  element  program. 

Instructions  for  SSCONDUCT 

SSCONDUCT  was  set  up  to  be  easily  modified  to  handle  new  problems.  Running  a  new  prob¬ 
lem  simply  involves  editing  subroutine  SSDAT A  by  following  instructions  provided  by  comment 
statements  in  the  program  to  adjust  the  variables  and  arrays.  The  variables  are  defined  at  the  be¬ 
ginning  of  SSDATA. 

However,  before  attempting  to  alter  the  data  subroutine,  the  user  should  make  a  drawing  of 
the  problem.  For  a  problem  with  more  than  one  thermal  conductivity,  a  number  should  be  assigned 
to  each  different  material  in  the  problem,  starting  with  1 .  Let  A  be  the  number  of  different  mater¬ 
ials  in  the  problem.  Assign  the  number  A  to  the  material  which  appears  the  most  in  the  problem. 

For  each  material,  obtain  the  thermal  conductivity  in  W/m  K.  If  a  material  is  located  on  a  con¬ 
vective  boundary,  obtain  the  convection  coefficient  in  W/m2  K. 

Now  the  finite  difference  grid  for  the  problem  should  be  drawn.  The  boundaries  of  the  grid 
should  match  with  those  in  the  problem  as  closely  as  possible.  Except  for  semi-infinite  bound¬ 
aries,  the  nodes  should  be  equally  spaced.  At  this  point,  a  decision  must  be  made  on  a  reasonable 
nodal  spacing.  To  date,  there  is  no  general  way  of  determining  the  optimum  nodal  spacing,  but 
there  are  the  following  guidelines: 

1.  The  nodes  should  represent  the  materials  in  the  grid  as  accurately  as  possible,  i.e.  no  node 
should  be  composed  partly  of  material  1  and  partly  of  material  2. 

2.  Regions  with  a  steep  temperature  gradient  are  best  represented  by  a  dense  grid. 

3.  In  general,  accuracy  decreases  as  internodal  spacing  increases. 

4.  The  more  nodes  used,  the  more  computer  storage  space  is  required,  and  the  more  computer 
time  will  be  needed  to  solve  the  problem. 

Bear  in  mind  that  the  use  of  semi-infinite  nodes  introduces  an  error  that  is  proportional  to 
the  area  of  the  semi-infinite  node.  Therefore,  the  semi-infinite  condition  should  be  used  only  in 
regions  of  small  temperature  gradients  when  knowledge  of  the  temperature  distribution  is  not 
critical. 

Let  y  be  the  number  of  rows  in  the  grid,  and  let  x  equal  the  number  of  columns  in  the  grid. 

The  grid  is  represented  in  SSCONDUCT  by  RAY  (I,  ),  K),  where  I  is  increased  from  1  toy,  J 
from  I  to  x,  and  K  from  I  to  3. 

It  may  be  desirable  to  make  three  drawings  of  the  grid.  The  first  time,  label  each  node  with 
temperatures  where  they  are  known.  This  will  represent  RAY  (I,  J,  1).  The  second  time,  label 
each  node  with  its  location  type,  as  defined  at  the  beginning  of  subroutine  SSDATA.  For  example, 
if  the  top  left-hand  corner  of  the  grid  is  at  constant  temperature,  label  it  "11”  (RAY  1 1, 1 , 2|  =11). 
The  grid  labeled  in  this  fashion  represents  RAY  (I,  ),  2).  These  labels  are  used  in  the  main  program 
to  assign  the  proper  equation  to  each  node.  The  third  time  the  grid  is  drawn,  label  each  node  with 
its  material  number,  as  determined  above.  This  will  represent  RAY  (I,  ),  3). 

A  copy  of  the  program,  including  SSDATA,  is  listed  in  Appendix  A.  Editing  SSDATA  to  agree 
with  the  new  problem  will  set  up  the  data  file  when  SSDATA  is  run. 
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When  locating  isotherms  in  the  final  temperature  distribution,  subroutine  ISOTHM  assumes  a 
uniform  internodal  spacing.  Thus  if  semi-infinite  boundaries  are  used,  edit  ISOTHM  as  directed 
in  the  comment  statements  at  the  beginning  of  the  subroutine  to  avoid  error  in  interpolations. 

When  SSCONDUCT  is  run,  it  asks  two  questions  of  the  user.  The  first  is  whether  or  not  sub¬ 
routine  SSDAT A  should  be  run  (it  should  be  run  if  the  formatted  data  file,  STSDAT,  is  empty  or 
if  the  user  desires  a  change  in  that  file).  The  second  question  is  whether  or  not  the  user  wants  sub¬ 
routine  ISOTHM  to  be  called  to  locate  isotherms  in  the  final  temperature  distribution.  File  STDOUT 
is  the  output  file  which  lists  the  initial  data  and  boundary  conditions  along  with  the  final  temper¬ 
ature  distribution. 

At  the  time  of  publication  of  this  report,  all  boundary  conditions  are  operational  on  all  sides 
of  the  grid,  except  for  the  semi-infinite  conditions,  which  are  operational  on  the  left  edge,  bottom, 
and  right  edge  of  the  grid  only. 

Instructions  for  FEHEAT 

The  first  step  in  using  program  FEHEAT  is  to  divide  the  region  to  be  modeled  into  triangular 
elements.  Care  must  be  taken  here  in  order  to  obtain  an  accurate  and  economical  solution.  Points 
which  will  help  are: 

1.  Use  the  smallest  elements  in  the  areas  where  you  would  expect  the  highest  temperature 
gradients. 

2.  Use  large  elements  where  small  gradients  are  suspected  in  order  to  minimize  the  number  of 
elements  and  minimize  computer  core  requirements  and  running  time  (i.e.  cost). 

3.  Avoid  the  use  of  triangles  with  high  aspect  ratios,  that  is,  triangles  with  sides  varying  greatly 
in  length. 

4.  When  several  materials  are  present,  approximate  their  boundaries  as  closely  as  possible  with 
the  element  boundaries.  Each  element  can  only  consist  of  one  material  in  the  model. 

The  next  step  is  to  assign  numbers  to  all  of  the  nodes.  Great  care  must  be  taken  here  to  ensure 
computational  efficiency.  The  form  of  the  matrix  of  equations  which  yields  the  solution  is  directly 
affected  by  the  numbering  of  the  nodes.  The  smaller  the  bandwidth  of  this  matrix,  the  lower  the 
computational  effort  required  to  solve  it.  To  minimize  the  bandwidth,  number  the  nodes  such  that, 
for  the  entire  region,  the  maximum  difference  between  any  two  node  numbers  of  any  element  is 
as  small  as  possible.  For  instance,  it's  better  to  have  a  maximum  difference  of  10  among  the  node 
numbers  of  several  elements  than  to  have  one  element  be  15  greater  than  the  smallest  element  and 
have  a  less  than  10  difference  among  the  remaining  elements.  The  element  or  elements  with  the 
largest  difference  between  node  numbers  determines  the  bandwidth  which  must  be  carried  for  the 
entire  solution. 

Once  the  nodes  have  been  numbered  the  elements  themselves  may  be  numbered.  They  may  be 
numbered  in  any  sequence  without  affecting  computation  efficiency.  Now  we  can  construct  the 
input  data  files.  They  are 

ED  =  clement  data  file 
NPD  =  nodal  point  data  file 

BCT  =  constant  temperature  boundary  conditions  file 

IHG  =  internal  heat  generation  file 

SHF  =  specified  heat  flux  file 

CBS  =  convective  boundary  segments  file 

QUAN  =  control  data  file 
TIOT  =  desired  isotherm  temperatures. 

Of  these  eight  input  data  files  only  ED,  NPD  and  QUAN  are  necessary  unless  the  boundary 
conditions  require  the  use  of  BCT,  IHG,  SHF  or  CBS  and  the  desire  for  isotherms  requires  the  use 
of  TIOT. 

The  input  data  for  ED,  the  element  data  file,  is  shown  for  a  sample  problem  in  Table  1.  The 
node  numbers  for  a  particular  element  must  be  entered  in  counterclockwise  order  around  the 
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element.  The  material  type  must  be  assigned  to  each  element  and  the  material  thermal  conductivity 
defined  in  the  program.  Any  consistent  set  of  units  may  be  used  for  the  thermal  conductivity, 
temperature,  and  position  variables.  The  numbers  in  this  file  must  be  in  the  format  (15,  416). 

NPD,  the  nodal  point  data  file,  is  constructed  exactly  as  shown  in  Table  2.  The  X  and  Y  coordi¬ 
nates  of  each  node  must  be  determined  from  the  drawing  of  the  region.  These  are  entered  in  the 
format  (15,  2F10.4). 

Constant  temperature  boundary  conditions  are  input  in  file  BCT.  The  format  of  this  file  is 
(boundary  condition  number,  node  number,  assigned  temperature).  These  may  be  listed  in  any 
order  of  node  numbers,  although  the  boundary  condition  numbers  must  be  consecutive.  A  typical 
BCT  file  might  appear  as: 


1 

3 

50.0000 

2 

7 

1 00.0000 

3 

10 

125.0000 

4 

2 

150.0000 

5 

39 

20.0000 

6 

8 

75.0000 

Thus,  for  instance,  boundary  condition  3  would  be  on  node  10  and  it  would  have  an  assigned 
temperature  of  125.  Again,  the  units  of  temperature  need  only  be  consistent  with  their  use  else¬ 
where.  Headings  on  the  printout  indicate  temperatures  arc  in  °F,  but  this  may  be  ignored  il  °C 
arc  used  for  input  throughout.  The  format  for  input  in  file  BCT  is  (15,  16,  F  10.4). 

F  ile  IHG  contains  the  boundary  conditions  for  internal  heat  generation.  The  format  of  this 
file  is  (boundary  condition  number,  element  number,  assigned  heat  generation  rate).  Similar  to 
file  BCT,  these  may  be  listed  in  any  order  of  element  numbers,  but  the  boundacy  condition  num¬ 
bers  must  be  consecutive.  A  typical  file  could  appear  as: 
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20.0000 
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10.0000 
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17.5000 
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22 

3.3333 
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16 

20.0000 

6 

4 

35.0000 

For  instance,  boundary  condition  5  would  be  on  element  16  where  the  rate  of  internal  heal 
generation  would  be  20.  Consistent  units  must  be  used  for  heat  generation  rate.  The  format  for 
the  input  in  file  IHG  is  identical  to  that  for  file  BCT;  it  is  (15,  16,  F  10.4). 

File  SHF  contains  the  boundary  conditions  of  specified  heat  flux.  An  ideally  insulated  boundary 
has  a  heat  flux  of  zero.  If  no  heat  flux  is  specified  for  a  node  on  a  physical  boundary  of  the  sys¬ 
tem,  it  becomes  a  zero  heat  flux  boundary,  provided,  of  course,  that  no  other  boundary  condition 
is  specified  there.  The  format  of  this  file  is  (boundary  condition  number,  node  number,  specif  ied 
heat  flux).  Again,  the  nodes  may  be  listed  in  any  order  but  the  boundary  condition  numbers 
must  be  consecutive.  Consistent  units  must  be  used.  A  typical  SHF  file  is: 

1  6  5.0000 

2  1 7  27.0000 

3  24  6.6667 

4  4  32. ‘s  000 

5  5  50.0000 
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6  87  9.5000 


For  example,  boundary  condition  2  would  be  on  node  17  where  the  specified  heat  flux  would 
be  27.  Again,  the  format  for  input  into  file  SHF  is  (15,  16,  F10.4). 

File  CBS  contains  boundary  conditions  for  segments  of  the  boundary  where  convection  occurs. 
A  boundary  segment  subject  to  convection  is  defined  by  the  nodes  at  its  end.  The  format  of  this 
file  is  (boundary  condition  number,  node  number  at  one  end,  node  number  at  the  other  end,  con¬ 
vective  heat  transfer  coefficient,  ambient  temperature).  The  boundary  segments  may  be  listed  in 
any  order  but  the  boundary  condition  number  must  be  consecutive.  Consistent  units  must  be  used 


for  the  convective  heat  coefficient  and 

1 

6 

7 

1.0000 

20.0000 

2 

27 

3 

1.2500 

25.0000 

3 

12 

20 

2.3750 

55.2750 

4 

21 

20 

1.0000 

20.0000 

In  this  case,  for  example,  convective  boundary  segment  3  would  have  node  1 2  at  one  end  and 
node  20  at  the  other.  The  convective  heat  transfer  coefficient  there  would  be  2.375  and  the 
ambient  temperature  55.275.  The  format  for  file  CBS  is  (315,  2F8.4). 

The  QUAN  file  contains  control  data  for  the  problem.  The  format  of  the  file  is  (NN,  NE,  MC, 
M|C,  MCI,  MC2,  NIT),  where 
NN  =  total  number  of  nodes 
NE  =  total  number  of  elements 

MC  =  total  number  of  fixed  temperature  boundary  conditions 
MJC  =  total  number  of  elements  with  internal  heal  generation 
MCI  =  total  number  of  nodes  with  specified  heat  flux 
MC2  =  total  number  of  boundary  segments  subject  to  convection 
NIT  =  number  of  isotherms  desired. 

A  typical  QUAN  file  is: 

104  166  26  0  0  0  6. 

Thus,  this  problem  would  have  104  nodes,  166  elements  and  26  constant  temperature  boundary 
conditions,  and  6  isotherms  would  be  desired.  The  format  of  this  file  is  (716). 

The  final  file  used  is  TIOT.  This  file  contains  the  temperature  of  the  isotherms  which  the  user 
wishes  to  locate.  The  format  is  (IT,,  IT2.  IT3,  IT4...,  ITn(T),  where  IT,  is  the  temperature  of  the  /th 
isotherm  desired,  I  <  /  <  NIT.  A  typical  TIOT  file  is: 

10.00 

30.00 

50.00 

70.00 

90.00 

100.00 

The  format  for  file  TIOT  is  (F8.2) 

With  the  necessary  files  ready,  program  FEHEAT  can  be  executed.  The  output  is  self  explan¬ 
atory,  giving  the  resultant  temperatures  at  each  node.  If  isotherms  arc  requested,  a  set  of  coordi¬ 
nate  pairs  which  will  form  the  isothermal  line  will  also  be  output.  Input  data  are  printed  out 
along  with  results  as  a  means  of  checking  to  see  that  all  data  are  piopcrly  read  in. 
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If  problems  occur  in  the  use  of  this  program  two  likely  areas  should  be  checked  first.  These  are: 

1.  All  dimension  and  common  statements  must  be  made  large  enough  for  the  problem.  Cur¬ 
rently,  the  dimensions  are  set  for  a  maximum  of  175  nodes,  300  elements  and  a  bandwidth  of  30. 

2.  Most  often,  mistakes  are  made  in  the  input  data,  particularly  in  files  ED  and  NPD.  A  very 
easy  way  to  check  these  files  is  to  write  a  short  plotting  program  which  would  read  first  the  coor¬ 
dinates  of  all  the  nodes  from  the  NPD  file  and  store  them.  Then  the  program  would  read  in  the 
nodes  for  each  of  the  elements  from  the  ED  file  and,  using  the  coordinates  of  these  nodes  as  ob¬ 
tained  from  NPD,  plot  each  element  individually  on  a  composite  plot.  After  all  the  elements  arc 
plotted  it  should  be  exactly  as  the  original  discretization  of  the  region  was.  If  any  stray  lines  are 
found,  element  sides  missing,  or  incorrect  boundaries  drawn,  the  problem  in  the  input  data  should 
be  easily  located. 


CONCLUSIONS 

Two  computer  programs  have  been  developed  to  solve  the  steady-state  heat  conduction  equation 
under  a  variety  of  boundary  conditions.  The  accuracy  is  the  same  for  each  when  modeling  problems 
with  rectangular  boundaries  and  no  semi-infinite  boundary  conditions  are  modeled. 

SSCONDUCT,  the  finite  difference  program,  is  by  far  the  e«s  lest  to  set  up  and  run  lor  a  new 
problem. 

FEHEAT,  the  finite  clement  program,  has  the  advantage  of  being  able  to  use  elements  of  varying 
size  throughout  the  problem.  This  provides  for  a  better  definition  of  curved  boundaries,  and  makes 
the  problem  cheaper  to  run  (less  computer  time)  if  large  elements  are  used  where  the  temperature 
gradient  is  small. 
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oooo  noon 


APPENDIX  A:  FINITE  DIFFERENCE  PROGRAM- DOCUMENTATION  AND 
SAMPLE  INPUT  AND  OUTPUT 

Program  SSCONDUCT 


SS CONDUCT 

THIS  FORTRAN  PROGRA*  SOLVES  FOR  STEADY  STATE  2-0  TEMPERATURE  DISTRIBUTION 
RESULTING  FROM  CONDUCTION  HEAT  TRANSF)K.  HOUNOARY  CONDITIONS  MAY  INCLUDE 
CONSTANT  TEMPERATURE  if  CONVECTIVE  5URFACFS.  SEMI-INFINITE*  AND  CONSTANT 
FLUX  DOJNJAR  IES.  THE  GMD  MAY  CONTAIN  MANY  DIFFERENT  CONDUCTIVITIES. 

OATA  FUR  SSCONDUCT  IS  u  A  T  H  E  R  t bY  sSDATA*  WHICH  PUTS  IT  INTO  FILE  STSOAT. 
SEf  SSDATA  FOR  AN  EXPLANATION  OF  VARIABLES. 

FOR  DIMENSIONS:  RMAT < Y«X ,2X*1 ) * F'VCl < Y* X, 1 > ,KXL ( Y* X. <X* 1 )) 

IMPLICIT  IiTEUE.iCA,:,.,  .U.X.Y,/) 

C  0  MR  0  N / M  1  I  /  A.X.Y.NISO 

COMMON /MIR/  t)S»  Jl  *T"AX*TMIN»F  I  <  A)  *  THK<  ?  »  ,H<2  )  ,  T  IS<H-  J  « 
f  RAY(1 DC. 130. 5  > 

COMMON /Mi/  K.MATUC0CA»c’0j),KVCT(10:uR.1'.  >»RXL(i;CCA,IOF) 

JI-'FNSU.N  IFLAGIAJ 

CALL  CONTRL <5 .‘STSU-T '»5> 

CALL  CONTRL  I •’.*STDO>iT,.G> 

«  R  I  T  L 11 •  1  ) 

1  FOn,,AT(lX.*IF  YOU  S  ITED  ThE  DATA  SUBROUTINE  SINCE  YOU  LAST  *. 

r  •  RAN  THIS,  TYPE  "l"*./»lX, • AND  ScDATA  WILL  BE  E  XE  C  UT  E  0  .  •  ,  /  .  1  X  , 

C  • D  f  H  L  R  « I  0  E  »  TYrE  ".  »•) 

K E  4  0  < 1 • »  «  E  rt  r : U  11  >  IT 

uo  to  : 

1111  CONI  IN Ji 

-  -1  1  T  r  (  1  ,  j) 

J  FORMAT  (  IX.  •lk.ROr  On  INPUT  FOR  YOUR  R  f  SPD’iSt  *  ) 

GO  TO 

2  CJNTINj. 

IFIIT.Ii.J  GO  10  a 
call  ssdata 

SO  1  U  u 

a  CO.TINJ' 

KLA.<5,lrJ>  „%»Ol 
10  FORMAT  (  lA.rP'1.".) 

cv  E  t  i  (  5  ,  ^  C  )  4  ,  X  ,  Y  .  N  I  .  0 

22  FOw f mT (  1  x,  a  1  j) 

UHK  CL  >.L  =  1  ,A> 
ix  E  A  f  (  5  ,  x  o  )  (  H  I  L  >  ,  L  =  1  ,  A  ! 

K  £  A  ."•  <  t> «  5  0  )  (  T  I  j  u(  ..  >  ♦  •.'=  1  »*J  I  :.C  > 

.5  C  FORMAT  Ilx.FH.Al 

kEmO<5,H)  (FMIMrI.O 
51  FO-x  ,-AT  <  1  A  ,  AF  <  .  A  ) 

ix  E  A  D<  b  ,  55  )  (  (  K  A  Y  (  I  ,  U,  5  )  ,  U=  1  ,  X  >  ,  i  =  1  ,  Y  ) 

REA. ( l , Mb)  <(RAY(I,J,D),U=1,X>,I=1,Y) 

R  L Au  <5 ,55)  (<nAY<I,J,l),U=l,X).lR]  ,Y> 

55  FORMAT < lx, 1  IF  7. 2) 

fa  CO  j T I f J U F 

WRITE  DATA  FOR  RUN  INTO  STUOUT  *»«••  *«»•  ... 

wRIMC<b,5CI 
W  R  I  T  C  <  1  *  b  0  > 

50  F  OK  M  A  T  f  1  X  »  •  OATA  I  OR  T  H I ' .  RUN  OF  SSCU*1  DUCT  :*  «/,/  ) 

WRIT  E(t>, 511  X 

WRI Tl< 1  .51  )  X 

51  FORMAT  < 1 X, »  X- *,  I  5  > 

WixlTE<6,52>  Y 
wRITE<1.52>  Y 

52  FORMAT! IX, »Y=*. 15) 

UR  I Tf  <6,5A>  A 
URITEll.bA)  A 

5 A  FORMAT < 1 X, »A=  •» 1 5  ) 

WRITE <6,56)  NISO 
WRITF<1,56)  NISO 

56  FOP  MAT < IX, »NI SO-' • I  '  ) 

WRITECw.B)  US 
WRlfE(l.H)  uS 

A  FORMAT < 1 X, *DS=» ,F 7.5) 
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MMOO  OOCJtJ 


58 
53 

59 
55 

57 


16 

15 


yRITE(t>,58)  0 1 
WRITEC1,58>  01 
FORMAT<1X,»OI=',F9.*> 
WRITE<6,53> 
yRITE!l»53) 

F0RMAT!1X,*TIS0!8>,6=1,NIS0:*> 
URITE<6»59)  !T1S0I8>«B=1,NIS0) 
URITE!1,59)  < TIS0!B)»B=1,NIS0 > 
FORMAT <1X,F7.2) 
yRITE<6»55) 
yRl TE ! 1 *55  > 

FORMAT! IX. »THK  <L> ,L  =  l, A: •> 
WRITE<6,30)  <THK(L),L=1«A) 

WRIT  £(1.30)  <T-tK<L),L=l»A> 

yRI  TE!6,57> 

y<UTE«1.57) 

format<ix,*h<l>,l=i.a:*) 
yKI  Tf..<6,30  >  (  HC  L)  »L  =  1  *  A) 

“  <  HIL  >  »L  =1  »  A  > 


y  R  I  T  E  Cl  .  3  0  > 
00  15  1=1.* 
y'RITE<6,16) 
WRITE!  1,16) 


<<*AY1I,U,3>,U=1,17),I 

<  <R AY!I .0,3) «U  =  1H,X> , I 

<  (w  AY  U,  J,2)  ,  0=1 , 17)  ,  I 
I(kAYII,0,2),0=18,X),I 

<  <  u A  Y (  1,0, 1) ,0=1 ,1 7) ,  I 


65 

b6 

67 

bR 

C  ** 


6* 


3  ,  F  1  (  I  > 

I  ,F  1(1  > 

F0RMAT!1X.*FI (,,11,,)=*,F9«*> 

CONT IN Jt 
WRITE(b»65) 

WH  I  TE! 6,990) 
yRlTE!b,892> 

•  RI  TCCo.893) 
yRITE<6,6b) 
wRITt!b,«93) 
yRI TE  C  6,692  > 

WRI TE  <6,893  ) 
yRITE(b,67) 

WR  I  Tt  <  6, 69  0 ) 

■RITE<6,892) 

■RITC<o.893)  <(r\AY<I,0,ll,0  =  lfi,X),I 

FORMAT  </, IX, *KAY<  1,0,3):*) 

FORMAT  </, IX, *KAl <  1,0,?) : ») 

FORMAT!/, lX.’hAYCI, 0,1):*) 
y  R  I  T  E  <6,bB) 

FORMAT  </, IX,* BOUND ARY  CONDITIONS:*  ) 
«*»»  *«* 

DI2=2.*DI-1. 
ld«=9«N0yIDT9 
IJy=2*x*l 
K  N  T  =  X  *  Y 
1)0  60  1=1, KNT 
oo  oi  l  =  l ,  I c y 
R  M  A  T  <  K  ,  L  )  =  j 
kVCT!K,1)=& 

CONTINUE 

CONTINUE 

00  b*  1=1,* 

IFLAG!K)=0. 

CONTINUE 


=  1,  Y) 
=  1,  Y) 
=  1  ,  Y  ) 
=  1  ,  Y  ) 
=  1  .  Y  > 
=  1  ,  Y  > 


C(7BP  T?!T7T2T'T5  FTJiTW  cnrFFs-ffF-rsiJ*TTo  NS  FOR  EACH  NODE 
FORM  MATRIX  RMAT  TO  CONTAIN  DIAGONAL  flTS  ONLY 
<00\T  =  f' 

00  2  01  I =  1  , Y 
00  2 '0  0=1 , X 

KOUNT  =  K  0  U  N  T  +1 
lNUtX=R AY! I ,0,5) 

1 F  < I .EQ. 1 )  GO  TO  11 
IND1=RAY< ( 1-1 >,U,3) 

TK1=2./!1./THK< IN01 Ml./THK!  INDEX)  ) 

11  CONTINUE 
IFtU.EO.l)  GO  TO  12 
I N02=KA Y< I , (0-1 > , 3) 

TK2=2./ U./THK! IND2>*1 ,/THK! INDEX)  ) 

12  CONTINUE 
IF(I.EQ.Y)  GO  TO  15 
I  NO  3  =  R  A  Y  < ( I ♦ 1 > ,0, 3) 

T*3=2./<  l./THK!  IN03  1*1  ./THKUN9EX)  ) 

13  CONTINUE 
IFIJ.EO.X)  GO  TO  1* 
l.ND4  =  R  A  Y  <  I  ,  <0*1  )  ♦  3) 

TK*=2./<  1  .  /  THK!  IINO*  »*1  ,/THKC  INDEX)  > 

1*  CONTINUE 

KKKrWAY! 1,0,2) 

X  I  =  X*I 
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X2=X*2 

XNM2*X)*1 

GO  TO  <  1  C  1  .102, 103,104, 105, 106, 107, 108,109, 110. Ill, 

4  1 12, 115*114,  US,  11  (.,117,118,119,  120, 121, 122, 123,124, 

4  12‘j,  1 2b,  1 2  7,  1 2t',  12s.  130. 131 . 1  32,  1  33)  ,  KKK 
C  INTERIOR  NODES 

101  CONTINUE 

RMAT (KOUNT  «  1  >  =  T  K 1 
R«AT<KOUNT.X>M ft? 

RMAT (KOUNT, XI >  =  -l . *  I TK 1 *TK 2*T K ? *T K 4 > 
ft  H  A  T  (KOUNT  »  X?  »  =  T  K 
RMAT  (KOUNT  ,  XN  >  =  TK3 
R VCT (KOUNT, 1 >  =  0. 

GO  TO  200 

C  CONSTANT  TEMP  UOJNOARYs 

102  COM  IN JC 

RMAT (KOUNT « XI  >  =  1. 

KVCT(KOUNT,l>=K«Y(I,J,I) 

2  GET  IFLAG  FOR  CONST  dOUND 

IF  (I.NE.U  00  TO  2i'2 
1FLAGI 1)=2 
60  TO  200 
202  CONTINUE 

JF(I.SE.Y)  bO  TO  ?r>, 

1FLAG( i) =2 
CO  TO  2UD 
23  3  CONTINUE 

IF  (  J.Nt  .  1  )  UJ  TO  22  <• 

IFlAG(2>=2 
GO  TO  200 
204  CONTINUE 

IFEAG(4)=? 

GO  TO  2  2  D 

C  CONNECT I  V  £  SJkFACE 
C  TCP 

134  cunt  1  n  u  t 

I  F  (  I  ,NE  .  1 )  C-0  TO  14 

I  F  L  A G  (  1  )-<* 

CONS  =  2  .  *H(  INDEX)  *'US 
ftKAT(K'JUNT»X>=TK2 

KNAT(K3UNT,Xl>  =  'l.*(TK2*2.*TK,i*TR4*"0N>) 

RMAT  (O’GNT  ,  *2  >  =  1<4 

RMAT  (KOUNT  ,  X'l  >  =  ..  .  *TK3 

R  V  C  T (KOUNT ,1)--1»LOnS»FAY( I, 0,1) 

GO  TO  2 ‘ G 
145  CONTINUE 
C  LEFT  SIDE 

I  A  (  J.Nt  .1  )  *>D  TO  I  «  ■ 

1 F  E  A  G  (  2  )  =  4 

C0VS  =  2.  *1>S*H(  INEEX) 

KMAT  (KOUNT ,  1  ) - T  K I 

*RAT  (KOUNT  ,  Xl  )--l  .  ‘  (  TM  »TK  1  .  «Tft4  ♦:  0.\  > 

R  N  A  T  (Kx>uNT*X?)  —  «  »  *  T  r.  4 
RMAT (KO JNT , XN  >  =  IK  ft 

RVCMKOjNT, 1>=-1. *C.  NS  *R  A Y ( I , U. 1 > 

GO  TO  203 
140  CONTINUE 
C  d  0  T  T  3  M 

IF(I.NE.Y)  GO  TO  141 
I  F  L  A  b  (  A  )  =  4 
C  0NS=2. *US  *H(  I  ME  x  ) 

RMa  T  (KOUNT  ,  1  )  =2.  *T».  1 
rx  '  A  T  (KOuNT  ,  X  )  -  T 
R  M  t  T  (KOUNT  ,  XT'  )  r  1  K  <1 

RMAT (KOuNT, XI  >  =  -l.*(2.»TKl*TK?,|K4,'0N’  ) 
RKC.T(K0UNT,1)=-1.,C^NS»RAY(1,J,1) 

GO  TO  2  3  u 
I  4 1  CONTINUE 

c  x  1  jm1  -j  ic r 

IFLa6(4>-4 

CONE =2, »0S  *  H (  I  \ c t  X  ) 
ft  MAT (KOUNT ,  1  ) :f < 1 
n MAT (K  GUNT  *X>=2.»Jr> 

ft  M  A  T  ( ft  3  U N  T  »  X 1  )  i  - 1  »  *  (  T  K  1  ♦  T  K  '  ♦  2  •  *  T  ft  ?  ♦  f.  ON  S  ) 

K  v  A  T  (KOuNT , X  N  )  - T  K  ■> 

KVCr(K&JNT,l»  =  'l.»f.N’>»RAV(I,J,X) 

GO  TO  ?■'  2 

0  CONSTANT  I T  E  N  I  Oft  i.OOf  - 

135  CONTINUE 

*RI Tl (6,151 >  I ,J«RAY(I • J,1 ) 

INI  FORMAT  (  1X,»C0NST  A  NT  INTERIOR  NODt  H  A  Y  (  *  ,  1 2.  ,  •  ,  •  ,  I  ?  ,  •  ,  •  ,  •  J  >  • 
i,  •b,,Fb.2) 

ft  14 A  I  (K  JUNT  ,  XI  )  :1  . 
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kvcu*ou*t.i»=R*¥*1*j,1> 

60  T9  f^.T  MUX  MOUNDAK* 

-  CONSTANT  HE  AT  FLU* 

'toi  continue 

'  IF  <  I  .ME-  •  1  >  CO  I  "  ^ 

iSkfcM^T.*»  =  TK2 

k.^aH10Unt.a1>-  i; 

jo  to  vrco 

lAt  CONTINUE 

. — 

kujix'uU.x* 

^AtTft<;.ONT*^»;TKH  >s 

T  ck^j  m  *  **  -  1 

v,0  TD  4.*'^ 

,  „ 

iff.U.Vr.l)  >' 

HnOUNI  ’  ..  .  <  Hi  >♦<  I<  ')♦<-*  '  7  "  '■*  *  ’ 

„■)  TO  i'  ■- 

\  (  *•  CONTINUE  f 

iCVi^o)f.T%i»  =  r^ 

,<  O ,  <  rt  r.  U  N  t  *  *  >  - ; :  * T/;  1 < T  r ! >  ♦ ,  , .  *  T '  u  >  *  ‘  '*■  ' ' 

k:';tu.iu<  t  : t<* 

^“nouNU^-ru.T.o, 


I  o 


'  7 


,T  1  N 7 


i-i-.'cr  'V 


r  J  J ' 


JA  -  T 


; ,  J 


i  > » *  t  t  < 


) 


i  jsi 


5  c  (  J  1  > 

TV 

4 ■jit-  n  *  .  > 

o.i  t  i  •  r-' 

(.  *1  *  Vv  T  1  *  •  r  I  1 

IF  <  J.Nt  .T>  •’°  " 

IFLA  oU  >|_j  . 

|  *  x  *  ^  -  1 1  I  »  «  T  <  >  V  '1  ‘  '1-  -!<1+  I-*1''- 
,  -U  T  <<  0  I  vT  »  xl  1  -  i 

UJ  ro  ?  • 

C0M1-.J'  ,  T  i  ■  ;• 

iflmI 

3 v„T  ( K J JNT »X\  1  -  x  • 

\  Hi  T  <H  VUN1  *  XT  >  =  '-  I  *  J  1  ^  r  j  ,  ,h;.y  (  U  J'  I  >  ’ 
RVC  T  (KJ'JNT  »1»  "  *  * 

c.0  to  2  '  'T 

.*•>»•'  «nRtt 

T  n  l  G  H  T 

li-«  =  THKlIM  tx)  Hl 

A  X  a  T  UOUNT  *  J  >  ^S,^/:)IOlMAl*l  1  r '  T 


ta 


O.TM/Ol  ) 


1 


KVCKKOUNT  *1>=-1.*(TKR/Dl  )»RAY(  I»J«1> 

GO  TO  2C0 

1000  FORMAT < 1 X* *KAY(  •,  12, • »*,  12,  • >  HAS  MO  EQUATION*) 

107  CONTINUE 

C  NODE  ADJACENT  TO  LEFT  SUE  SEMI-1VF  UOUNOARY 

RMA  I  (KOUNT  ♦  1 >  =TK  1 
RMAT (KOUNT  , X>  =TK2/Q1 

RMAT (KOUNT  *  X 1  )-“l •*( TK1*TK?/0I*TK3*TK4) 

RMAT (KOUNT  *  X2  >  =  TK4 
*.MAT  (KOUNT  *  XN  )  -  T  <  M 
«VCT(KUNT*1>=0. 

GO  TO  ?rC 
103  CONTINUE 

C  NODE  ADJACENT  TO  BOTTOM  SE^I-INF  UN  OR  V 

RMAT ( K 0 JN T  «  1) =  TK1 
KMA1(K0U\T»X)=TK? 

K«AT(K0Jf4T*Xl)  =  -l.*(TKl*T<r*TK5/DI*TK<t) 

R MAI  (KOUNT  »X?  >-TKA 
K  M  A  T (KUuNT . XN )  =  TK  5 / '  <  I 
RVC  I  (K  Tjr.T  «  1  )  =0. 

GO  TO  ?oo 
103  CONTINUE 

C  NODE  ADJACENT  TO  RIGHT  St M I - INF  IN  I  T E  RONDARY 

T  (ROUNT  ♦  1  )  =T.U 
RMAT (KOUNT *X)=IK? 

RMAT(K0UNT,X1)=-1.*(TK1*TKL*TKX*-TKA/0I) 

RMAT  (K0UNT*XL')-T;<4/  'I 
K'4-T  (KOUNT  iXN  )  -  t  k 
rv.:t<ount«1)  =  3. 

g  ■;  t  o  o :  o 

1 1  J  CONTINUE 

r  ,  u  l<  A-Tjtcr*  T  to  TOT*  Sf  I'.I  -  INF  IN  IT  E  BOUNDARY 

r  v r  < < o unt  « l  >  -  T r  l  /  j  i 

v  ••  I  (  K  j  o  N  T  *  X  )  -  r  \  ,  • 

-<  /,  A  T  (  K  'j  J  ,  T  |«1  1  -  -1  .  *  t  Tr]  /  r.  I  ♦  T  K  i  ♦Tr5*Tra> 

A  '  A  T  (  R  U  *  T  t  A.’  )  -  f  R  A 
R-.AT(KO.:jT,XN)-t<- 
h  (uT  (*  ,ij  iT  *1  )  U  • 
c  0  r  o  r  *  o 
in  ; r  l  v  i 1 

TON  *  I  ANT  ;  N*  - 
1  T  (R  ,'V'T  ,  »}  )  --  i  . 

R  .  T  <  R  0  J  N  I  ,1)-j*.Y(I  i  J  %  1  ) 

,  -  I  r  x  ( •  <  1 7 1 )  !ijf.;ni*j.D 

171  F  ;  V  •  il  (  l  X.  ♦  .  ON  -T  «\T  CORNER  A  Y  (•*  1  I  <»*,»,*  I  >  =  •  , 

>  r  n  •  „■  ) 

r  j  r  u  r  ■  J 
i  1 2  ;  ■  1 1  m  j  ? 

i  f  «  j  .  n  •" .  i )  m  t  j  ii  : 

IF(I.Nt.l)  JN  Ik  II  l 
C  TO?  U  U  C  r. C'.'  T  •  LJ  ' 

k  ■'  A  T  (  K  v  U  N  T  *  A  V  )  -  I  K  H 
RM  At  ( K  0  J  N I  ,  X  N  )  —  T  <  ' 

N  *  A  I (K  0 JNT  *  a  l  >  :-1.»(TM*T'h) 
n  V  t  T  (Kf.  jNT  ,  I  )  -  -  :  G  *  (  :  1  (  1  )  *F  l(  C  >  > 

.  u  i  t  f  ( t  , ;  i  7  > 

J  J  7  FOR  VAT  (IX*  *10?  LEFT  COR  Ml  K  DON'T  FlrC) 
c-o  To  LCj 
l l? I  CONTINUE 

:  •G1TCM  Li.FT  COtNfh  COMM  ?  I  JX 

«  HA T (K  DUVT  *  X?  >  =  I  K A 
RMAT  (KOUNT,  1  )  =  T  1 
•*"  *  1  (K  >U  NT  »Xl  )--l.*(TK‘-*TKi  ) 

«V>!(KOJNT,l):-rjt(l  1C  )«FI(  ‘  )  ) 

. •<  i  rr<F,,2;f  i 

UrH  FORMAT  (1X«*  '.'0  T I  iM  L  FT  CORNER  (  >-;ST  FL1'/*) 

u T  TO  2  J  0 
112'.  CONTINJ' 

I  f  (  J  ,!jr  .  T  )  r-c  T'i  i  I.  2 
C  'OTTOM  RIGHT  CG-xNFr  C  0  T  Ku* 

R  »  I  (MOUNT  *  1  )  -  T  M  1 
RMAT  (KOUNT, X)  :TK? 
kv“'T  (kuU\T»x1  )--1«*(T“I*Tk  ) 

R  YET  (MOUNT  *  1  >  =-.'S»  l  F  I  (  >  *F  I  (  A  )  ) 

«i  R  I  ft  (jkJR) 

200  F  jHMAT(  1X**:‘0IT  xM  RIGHT  C0*NtF  CwV.T  IL'JC) 

GO  TC  2  "  C 
112,.  CONTINUE 

C  TOP  RIGHT  C  HU  CONST  FLUX 

RMAT (KOUNT ,X)rT»: 

•(.MAT  (K  OUST  ,  XN  )  =  TK  A 
KM\T(K0UNT,X1)--1.»(TKU*-TK  > 

K  VCT  (KOUNT  *  1  )  =-  uS*U  I  (  1  >  *F  I  (A  )  ) 

K  R  I  T  f  (  h  *  2 1  0  ) 
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210  FORMAT ( 1 X  *  *  TOP  RIGHT  CORNER  CONST  FLUX»> 

GO  TO  200 

113  CONTINUE 

IF(U.NE.l)  GO  TO  1130 
IFU.NE.l)  GO  TO  1131 
C  TOP  LEFT  CORNER  CONVECTIVE 

C0NS=2.*DS*H< INDEX) 

RMAT (KOUNT  »  XI )=-l,*(TK3,TK4*C0NS) 

KMAT (KOUNT  *  XN >  =  TK  3 
RMAT(K3UNT»X2)=TK4 

R  VC  T (KOUNT ♦ 1 >  =  -l • *CONS*RAY (I »  J  *  1 ) 

WRITE(b»211) 

211  FORMAT ( IX* *  TOP  LEFT  CORNER  CONVECTIVE*) 

GO  TO  2 0 C 

1131  CONTINUE 

C  BOTTOM  LEFT  CORNER  CONVECTIVE 

C0NS=2.*DS*H{ INDEX) 

RM4  T (KOUNT  « 1 >  =TK  1 

RMAT (KOUNT .Xl >=-l .* <TK1*TK 4, CONS) 

RMAT (KOUNT  *  X2 ) -TK  4 

RVCT(KOUNT»l)=-l.  »C.jNS*K  AY  (  I  ,  U,  1) 
kRITt.  (ft*212) 

212  FORMAT (IX, ‘BOTTOM  LEFT  CORNER  CONVECTIVE*) 

go  to  zro 

1130  CONTINUE 

IF(I.NE.Y)  GO  TO  1132 
C  BOTTOM  RIGHT  CORNER  CONVECTIVE 

cong=2 .*os*h( index) 

R  MAI  (KC'JNI  «  1  )  =TK  1 
RMAT (KOUNT  «  X) = T K 2 

RMAT  (KOUNT*  XI  )  =  -l.*(TKl*TK?*CONr>) 

K VCT (KOUNT , 1 ) =-l . *U  NS*R AY ( 1 , J«l> 

W  r(  1  T  E  (  ft  »  2  1  3  ) 

213  FOR "AT ( 1X»*F0TT  0*  Kl&HT  CORNER  CONVECTIVE*) 

GO  1 0  2  'Jo 

1132  CONTINUE 

C  TUP  RIGHT  COaNER  CONVECTIVE 

C0NS  =  2.  *;JS*H(  I'.'JEX) 

R  v  A  f (KOUNT  ,  X) =T*i 

RMAT (KOUNT  , XI >:-l .» ( T K ? *TK3 ♦CONG ) 

RMAT (KOUNT  ,  XN  )  - TK  5 

RVC T (KOUNT ,1>=-1.*C0NS*RAY CI» J«l) 

m  R  I  T  E  (  ft  ,  2  1  4  ) 

214  F0RMAT(1X«*T0F  tUGHT  COR  NE  F  CONVECTIVE*) 

GO  TO  ?  0  C 

lift  CONTINUE 

C  CORNER:  VERT  =  CONVECT,HO«U  =  CONST  FLUX 

IF(I.EJ.I)  FH  1  =  F I ( 1 ) 

IF(1.E«.Y)  F  h I  =  f  ■  I ( 3 ) 

GO  TO  1 1G0 
117  CONTINUE 

C  CORNER:  HOR I2-C3NVECT,VtRT=C0NST  FLUX 

lE(J.EO.l)  PHUFIC) 

IMJ.LJ.K1  PH  I  =  r  1  (  4  ) 
lift?  CONI INJE 

IF(J.NE.l)  GO  TO  Uhl 
IF(l.NE.l)  CO  TO  llf.2 
0  TOP  LEFT  CORNER  CONVECT,  CFLUX 

CO  U  =  0$  *H( INDEX) 

RMAT (KOUNT ♦ XN  >  =  1<  3 
RMA I (KOUNT , X2  >:TK4 

R  M  A  T ( KOUNT , X 1  )  =  - 1 . » ( T  K  3* T K  4  *C  ONS  > 

K  V  C  T  (  K  0  U  N  T  ,  1  )  =  - 1 «  •  t ')  N  S  *  R  A  Y  (  I  ,  J  »  1 )  •  PH  I  *  D  S 
AHlTt(o,215) 

215  FORMAT  (  IX,  *TOP  LEFT  CORNER  CONVECTIVE  t,  0  FLUX*) 

GO  TO  2 " 0 

1 1ft 2  CONTINUE 

C  dOTTOM  LlFT  CORNER  CONVECT,  I;  FLUX 

CONS  =  DS*HC I NiU  X ) 

RMAT (KOUNT , 1 ) =  T  K 1 
RMAT (KOUNT  »X2  )  =  TK4 

KMAT (KOUNT, XI  )  =  -l .» (TK1 ♦TK4*C0NF) 

RVCT  (KOUNT, l)=-l.»Ci'NS*RAY(I,U,l)-FHI»DS 
■RITE (ft, 21ft) 

21ft  FORMAT  (  IX,* BOTTOM  LhFT  CORNER  CONVECTIVE  &  CONST  FLUX*) 
GO  TO  2  '1 0 
llbl  CONTINUE 

IF(  I.NE.Y)  GO  TO  lit  3 

C  BOTTOM.  RIGHT  CORNER  CONVECT,  C  FLUX 

CONS  =  OS*H(INDtX ) 

RMAT (KOUNT , 1 ) =TK1 
KMAT (KOUNT ,X  >  =T  K2 

RMAT  (  KOUNT,  XI  >  =  -l.»(TKl*TK.'*CONS) 

R VCT (KOUNT, 1 )=-l.«CONS*PAY(I,U,l> 
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WRITE(6»217> 

217  FORMAT  (  IX, ‘BOTTOM  RIGHT  CORNER  CONVICT  <,  G  FLUX*) 

GO  TO  200 

1HJ  CONTINUt 

C  TOP  RIGHT  CORNER  ECNVECT.C  UJX 

CONG=DS*H( l NOE  X ) 

RMAT (KOUNT ,X>=TK2 
RMAT (KOUNT  ,  XN  >=TK  3 

RMAT (KOUNT, XI >  =  -l .  * (TK2,TK3*C0NS1 

RVCT (KOUNT , 1 >=-l ,»C  ONO** AY ( 1 , U, 1 >-PHl*rc 

KRlTE(G*21b> 

218  FOR  MAT <  1  X, ‘TOP  RIGHT  CORNER  CONVECT,  CONST  FLJX»> 

GO  TO  200 

11^  CONTINUE 

IF(j.NE.l)  GO  TO  1140 
IFd.NE.  1)  GO  TO  1 1 A  1 

C  TOP  LEFT  COrNEF  Si.  “I-INTINITF  ON  20TH  SIOFG 

W»ITEd,lPCC>  1,0 
GO  TO  79 8 

1141  CONTINUE 

:  BOTTOM  LIFT  COftNTR  SEMI-INFIfdTF  ON  BOTH  GIGE 

K “ A T ( MOUNT , 1 ) =T  v  1 
RMAT (KOUNT *  X  2  )  =  1  X  4 
T  K  1  -  T  X  1 
T  X  L  =  T  X  1 
T  X .  -  T  X  1 

TLF T=RAY( I , J.l) 

To  I "  =  R  A Y  < I , J,  1  > 

01=100. 

HUM  (KOUNT  ,  1  )  =  T  “  1 
KHAT(KJONT,X=)=TX4 

K-iA  T(X0ONT,xl):-l.*(TXl+TX4,Ti<L4TXr> 
rVC  T (KOUNT ,1)=-TKL*  TLF  T-Tx  *  T 1  T  M 

■  rite  <  b » 2  r ;  > 

>21  FORMAT (lx* ‘BOTTOM  L_FT  CORNER  JZ  «  I  -  I  NF  I  N  l  T  T  •  ) 

CO  TO  2^0 

1142  CONTINUt 
IFd.NF.Y)  00  To  1142 

2  .-UTTOX  KhHT  CORNlR  2  E  M I -INF  { .4  1  T  ON  ••  0  I H  ST::'S 

rm.-.  T  (KOUNT  tlUTM 
R“hI CK0UNT,X)=TK2 
TK~=  TKl 
TKa=TK1 

TR I T=RA Y( I , J,l> 

T,)l  "-R  A  Y  (  I  ,  J,  1 1 
01=100. 

RM  f  T  (KOUNT  ,  XI  )  =  -l  .  *  (T*  1,TKJ,TKK*TK  •  » 

R  V  C  1  (  X  0  U  N  T  ,  1  >  =  -  T  K  R  «  t  R  I  f  -  T  K  *  T  B  T  w 
A  <\  I  T  f  (  6 , 2  2V  ) 

222  F'lRMATdX, ‘BOTTOM  RIGHT  COi'NEF  GE  M  I  -  1  NT  l  1  T  E  •  > 
f> To  2  2  0 
1142  CONTINUE 

C  TOP  RIGHT  CC.RNfR  f  1 -I  NF  I N I  T  E  ON  .OTH  SIfiS 

xRITF  (1,11100)  1,0 

fj  TO  .4.' 

11^  continue 

C  ;>  <  J  A  R  E  INTERIOR  N02E  ‘.‘EXT  TO  T  JO  St  KJ  -INF  SIOES 

IF(O.NC.2>  GO  TC  11  0 
I F (  1  ,NE  .2)  GO  TO  11:1 
C  TOP  LEFT 

A  x  I  l  E  (  1  ,  1  T  C  0  )  1,0 
GO  TO  '’c9 
1151  CONTINUE 
C  BOTTO*  LEFT 

«  X  A  T (XOONT , 1 ) =TK 1 
RMAT (KOUNT *X)=TK2/Ol 

KMAT(K0UNT,Xl)=-l.*(TKl*TKE/ni4T‘>3/r'I*Tl'.  4) 

RMAT (X JUNT  ,  X?  >  =  TK  4 
RMAT  (XOONT, XN  >  =  TK3/:  I 
RVCT (KOUNT,  I ) =0. 

C  JR  I TC(fc,22b>  I,J 

22 S  FORMAT ( 1  X  »  •  H  A  Y  ( ,»I2,,,,,I?,,2)  =  1N*  > 

GO  TO  200 
1150  CONTINUE 
YY=Y-1 

IF (I .NE. YY)  oO  TO  115? 

C  bottom  RIuHT 

RMAT (KOUNT, 1)=TK1 
RMAT (KOUNT »  X  )  =  7  !\  2 

R'1AT(K0UNT,X1»  =  -1.*<TKI»TK^,TK3/0I*1K4/0I) 

R  M  A  I (KOUNT  *X2 >=TK4/ul 
RMAT  (KOUNT  ,  XN)  =  TK5/'jJ 
RVCT (KOUNT  ,  1  )  =  G. 

GO  TO  200 
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1152  CONTINUE. 

TOP  RIGHT 
whlTE(l.lOOO)  I.J 
GO  TO  098 
CONTINUE 

BOTTOM  RIGHT  CORNER  HORIZ  S l M 1 - I NNF , VER T  CONST  FLX 
012=2. *01-1. 

T  K  0  =  T  K  1 

RMAT  (KOUNT  »1)=TK1/10I*2.) 

RMAT (KOUNT, X>=Tk2*D12 

RMAT (KOUNT  »X1 >=-l .*<  TK1/(D 1*2. > ♦TK?*DI2*TKB/(DI*2.»  > 
RVCT(K0UNT»1>=-1.*(TK9/([)I*2.>>*RAY(I»J»1>-FI(4>*DI2*0S 

GO  TO  2  C  0 
CONT INUC 

BOTTOM  RIGHT  SCUAKF  NEXT  TO  SL^I-INF  »RT  SIDC  CONST  F LX 
l)I2  =  2.*DI-l. 

RMAT (KOJNT  » 1 >  =TKl /2  . 

KMAT(KGUNT,X)=TK2 

RMAT(KOUNT,X1)=-1.*<TK1/2.*TK2*TK3/01«2.») 

RMAT (KOUNT  «XN )  =  T< 3/(01 *2.) 
rVCT  (KOUNT  .1  )  =  -E  1(41*0*. 

GO  TO  2C0 
CONTINUE 

BOTTOM  LEFT  SGUARF  NET  TO  St M ! -1 Nc  ,LEF T  SIDE  CONST  FLX 
U  I  •  »U  I  -  1  • 

RMAT (KOUNT , 1>=TK1 /?. 

R1AT (KOUNT  »X1 )  =  - 1 . * (TK1/2.*Tk4*TK3/(DI*2.)  1 
RMAT  (KOUNT  ,  X2  »  =  TK<. 

RMAT (KOUNT  «XN  )  =  T<3/(DI  *2.) 

RVCT (KOUNT ,  1  >  =-F  1  (  2  >  •  OF 

GO  TO  2 1  3 
CON  T  INJE 

-OTTOM  L  F  FT  CORNER  H  :K  I  Z  SE M I  -  1  NR . V F k T  CONST  FLX 
012  =  2. *01 -1  . 

T  < !'  =  T  K  l 

RMAT (KJUNT . I)=TKI/(I*2.) 

RMAT(K0uNT»X1>  =  -1.*(TK1/(OI*2.>*TK4*1)I2*TkU/(0T*?.>> 

K  MAT (KOUNT  «X2)  =  TK4»lI2 

R VCT ( KOUNT, 1)=-1. *(TK«/(D1 *?. ) ) *RAY(  1. J. 11 

GO  TO  2  v  C 
CONTINUE 

SEMI -INF  NODE  ABOVE  SEMI-InF  CORNER  NODE 
IF(j«Nr.i)  on  t;.  13'  o 
LEFT  olti 

TKJ  =  '»l  /  (  (  (  D  I -.5  >/THK  (  I  NO  3)  >  ♦(  .N/THK  (  INDt  X)  )  ) 
TK4=^I/({(DI-.r)/THK(IN0EX)>*(.5/THK(IN)4>>) 

RMAT (KOUNT , 1 ) =TK 1 *0  12 
rt  MAT  <  K  C  U  N  T  »  X  2  )  =  T  K  4  /  '  I 
RM4T(K0UNT»XN)=T< 3 *(012/01) 

RMAT(KOUNT.Xl>  =  -i.*(TKl*OI<r*TK4/OI*TK3*OI2/OI*THK(INOEX>/OI) 
RVCT(KOUNT»l)=“l.*(THK(iNOFX)/OI)*RAY(I.J,l> 

GO  T  U  2E0 

1300  CONTINUE 

RIGHT  SIDE 

TKL  =  01  /  (  (  (  J  I  -  .5)  /TnK(  IMTEX  »  )♦  (  ,-N/THK  (  I  \  r>2)  >  ) 

Tko  =  OI/( ( (C I-. ;)/Tn* ( I  NO 3) >♦ ( .S/THK(  IN~CX) 1 ) 

RM(  T  (KOUNT,  1)  =  TM  -  J  IP 

R«-T  (KOUNT  ,X)=Tk2/i)I 

R'MAT  (KJUNT  , XN  )  =  TK3*(UI2/0I  J 

KMA  T  (KOUNT  ,  XI  )  =  -l  .  *  (TK1*0I  ;:*TK.  /OI  *TKJ*N  I2/0I*ThK(  INOFX  )  /O  I  ) 
RVCT  (  KOUNT  ,1  >  =  -l  .»  *  THK  (  INDEX  )  /Ft  I  >  *R  AY(  l  ,  J.  II 

GC  TO  200 
153  CONTINUE 

SEMI-INF  NONE  TO  RIGHT  OF  SFMI-INF  CORNER  NODF 
IF(I.NE.l)  GO  TO  133(1 

TOP 

WKI TE(  1  ,1C?0)  I  ,U 

GO  TO  '">>< 

1330  CONTINUE 

BOTTOM 

T  K  1  =  C I  /  (  (  (  I  -  .  J  1  /  THK  (  I  NOE  X  »  »  ♦  (  .  f  / T  HR  (  I  N  J 1  »  )  ) 

TK2  =  D1  /  (  (  (  0  I  -  .N  >  /  T  HK  (  IN'UZ  >  >♦<  .‘j/THK(  IN~EX>  >  ) 

R  M  A  T  (KOUNT,  1  )  =  T  R  1  /(.  I 
RMAT (KOJNT ,X)=TK2»UI?/Ul 

RMAT  (KOUNT.Xl  )  =  -I.«(TKl/0I*TK2»0I2/0I*TK4*r)I2*THK(  INOEXl/OI) 
RMAT (KOUNT , X?  >  =  T K 4 • C I  2 

RVCT(K0JNT.1)=-1.*(THK(IN0EX>/UI»*KAY(1,J,1J 
GO  TO  200 

120  CONTINUE 

121  CONTINUE 

124  CONT  INUF 

125  CONTINUE 
12b  CONTINUE 
127  CONTINUE 


C 

118 

C 

11B 

-  123 
C 

^ 122 

130 
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129  CONTINUE 
129  CONTINUE 
1  i2  COM  I  NO  C 
131  CONTINUE 

WRITE(1»1303)  I.J 
GO  TO  996 

200  CONTINJt 

201  CONTINUE 


LOOP  TO  WRITE  iiOUNJARY  TYPES  ONTO  STOfu/T 
DO  r'G  K  =  l,  9 

IF  I  K  .  N  F  .  1  )  G  f  I  TO  9  0  1 
W  K  I  T  £  <  F.  »  9  C  0  9  ) 

900  9  FORMAT  I  1 X, *  TOP  UOJNjARY*) 

GO  T  li  ')  j  0  9 
9001  CONTINUE 

IF  (  *  .  N  L  .  :> )  GO  TO  9.02 
WKlTLlp,9GOb> 

9  0  0  (  F  OF,  I'  A  T  (  l  X  ,  *  Lt  FT  tlOuNGAr  Y  *  ) 

G  J  I  Ll  *  ■  l  4 

TOO r  CO.IINUI 

IF  (K.VI.i)  b-1  IG  -  2  3 

•  n  I  T  E  <  F  ,  9  0  'j  1 ) 

-037  F,,v"  AT  (  1  a,  *r,  )  T  i  ■  "  ■  .  J  \  i..  A  H  Y  *  > 

09  to  -.19 
1  0  0  i  L  G  E1  T  I  N  U  r 

*  III  (  oi  '  j  ,  ,1  i 

9  3  :  M  F  OP  •'  A  T  I  1  X  ,  •  9  1  Gr  !  .  ;  N  J  A  9  Y  «  ) 

:•  :■  0  cost  I  r*  u  e 

1 F I  I  =  L Am  »<  >  .NL  ..  >  ,  T  r  -l 

»  I  T  t  (  .  ,  -  9  ) 

9 I  9  ’ 


•1 


?  ' 


'  A 


"j 

9  7 

9  H 
1  6 


C'.'  T  I  Mi;' 

IF  <  I  F  L  A  G  I  X  >  .  ■«  i  .  •*  )  .9  TO  < 

*  \  I  T  L  I £  ,  -  7  ) 

GO  T  J  • 

1 1  •  i  I  N  J  , 

IF  (  ifl  >.g(<  >  .v„ .0  ro 

,  I  f  ’  (.»“  5  ) 

G'l  Tu  ••  ? 

CO  ,  I  I  Vu‘ 

IFtlrL*  ><<).’«►»  )  ,  T9  ; 2 

»  1  I  I  E  t  ,  9  b  ) 

C  J  TO  •  j 
CG9  I  I  a  J 
F  0  i  9  A  T  (  1  h  ♦  ,  • 

'AT  <  In*,* 

F  O  \  4  T  (  1 1 1  ♦  «  * 
f.Ga 

c  .9 1  ii  n  j 


(  9  x.  j  T  A  \  I  Tf.  VE-AIJO' 
L'WeCTIVC  >l)RF  ACE  *  > 
CONSTANT  PEAT  FLUX') 

S  r  *•  i  -  J  -,l  I  M  I  E  •  ) 


•  > 


X  L  11119*70)  I  •  A 

■jx  nr ti*7?)  it. 

7  2  F0-i'*AT</*lX.‘!  A'  -  Wl  .) T  H  =  •  .  I  ) 

:  S3  L  V  t.  -i  T  R  I  X  »"AT 
NIT  -x 
N  b  0  =  X 
N  -  Y  *  X 
I A  =  Y  »X 
«-  1 


I;-=Y*x 

i  jr> o 

I  E 9  - d7 9 

WR 1 T  E  « 1 *63  ) 

<«IU(tixO) 

90  FORMATIlX,*  N  L  G  »  NUC,  N *  IF*  ",  J-,JUiVi,  ILL 
VRITE<1»81>  NLC*NUC*N*IA,M»lt'*lJ3u*lLi> 

*Rl TE<b,«l 7  NLC»NUC,N,I A ,M,I9 , I JOH  *  I  IF 
31  Fort  MAT  I 1  X  »  8  1 5  > 


•  ) 


CALL  BANDMX«N*NLC*YJC*IA*M,I<*IJ3T,I{S) 
wRI  T  c.  I  1 ,7A  )  IErt 
WRITE (G, 79)  I  tR 

79  FORMAT!/, IX, • If p=*» 13, •: *,/,S x, *JF  IfR=l?  ,  MATRIX  IS  SINGULAR*, 

S  /,9X,*IF  IER=3,  EVERYTHING  IS  GRAY.*) 

IF  I IER .E J.124  >  GO  Tn  998 
K  =  J 

00  310  1  =  1, Y 
00  300  3=1,  X 
K  =K  ♦  1 

K Ar< I, J, 1>  =RVCT  IK  ,1  ) 

300  CONTINUE 
3)0  CONTINUE 

WRITE<8,891> 
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R91  FORMAT < IX, *FI ML  TEMPERATURE  DISTRIBUTION:*) 

RRITEIb.flGG)  CIRAYI  WJ,1),J=1,1  7),  I=1,Y) 
rfR  I  TE  ( f>  ,  b92  ) 

WR  I  Tt.  I  b  .893  )  <(KAY(I«J»1),J=18,X)«I=1*Y) 

890  F0RMAT(1X»17F7.2) 

392  FORMAT </.l X, *. .THE  REST  OF  THE  COLUMNS:*) 

393  FORMAT  < 1X.4F7.2) 

CALL  ISOTHM 

GO  TO  99fl 
993  CONTINUE 

jRITEU.9999) 

9999  FORMAT  IlX.’EHHDR  IN  ASSIGNING  NODAL  LOCATIONS*) 

99S  CONTINUE 

CALL  CONTPL I4.*STSDAT • .5) 

CALL  C0NTRL<4,»STD0UT*,6) 

CALL  EXIT 
END 

-r?rTTRff"FT^rT5~T5!5TPriTTisTnfr-55rcnr(Trr - 

SUBROUTINE  ISOTHM 

IMPLICIT  INTEGER  I  A.I-  .u.b.X.Y.Z) 

common/mii/  a.x.y.niso 

C  0  fi  M  u  N  /  M  I R  /  LS«0I»TMAX*TMIn*FI(4)*TtJK(?)*H<?)»TIS0<f5)» 
i  RAYC120.10C.3) 

DIMENSION  EX«Yt9.bO).  EYftYt9»bO>»  COUNT (9) 

01  Ml  NS  I  ON  IX(S)«IY<3) 

CALL  CONTRL <?.*POINTS* .7) 

CNTPTr? 

UO  3  L -  1  »  9 
C'JJNT(L)-0 
i  CONTINUf 

CHANGE  THE  NXT  T  »  D  ST^TS  TO  AGREE  mITh  DIMENSION: 

DO  u  K -  1  *40 
0  0  7  8  =  1,9 
EXrxY<H,K)  =  : 

£  YRY(u.-C)  =  E 
7  CO-JTINJL 

9  CONTINUE 

TMAX-TEMP  OF  HOTTEST  ISOTHERM 
TMIN-TEMP  OF  COLDEST  ISOTHERM 
NISO-  NO  OF  ISOTHERMS 

LET  1=  HOTTrST  ISOTHERM 

COUNT  <rl=COUNTfR  FOR  NO.  f.LTG  IN  EACH  ISOTHERM 
TISOIU):  TEMP  OF  lACH  ISOTHERM 
EXKY<H»<)=ARMAY  FUR  x -C OCHDI N AT [ S 
r.  INDICATES  «HlCH  ISOTHERM 

m=C0UNT(S)  a  INDICATES  POSITION  OF  PT  IN  ISOTHER**  LIST 
LOOP  TO  LOCATE  ISOTHERMS 

FOr  LEFT  SEMI -1  \NF  .LET  XL  =  ‘, 

FOR  -OTTO.M  SEMI-INF  .LE  T  Yh-Y-3 

FOR  RIGHT  Sl MI-INF  .LET  XR=X-3.  AND  COMMENT  OUT  LOOP  Pr 
FOR  TOP  SLMI-TNF.LET  YT  =  i,  AND  COMMENT  OUT  LOOP  BfeO 
IE  NOT  US  I N  o  G£MI -INF,xL-1.xR=X-1,YTp2,YR=Y 

XL=  i 
X  ?.  =  x  - 1 
Yb  = Y-3 
Y  T 

00  I'O  I  :  Y  T  ,  Y  H 
DO  ?jC  J  =  X  L  ,  X  r 

EXAMINE  TEMPS  HOR  I /'TNTALLY 
K J  =  R AT ( I , J»  1) 

R JI=RAY< I. < J* 1 ) . 1 ) 

IF C (KJ.GT . T ISO* 1 > > . AND. ( R J1 .GT . 1  ISO <  1  )  )  )  GO  TO  900 
I F < <RJ .LT .T ISOINi SO) ) . AND. CR J1 .LT.TI SOCNlSO > > >  GO  TO  5CC 
DO  SCO  3=1. NISO 

IFU  T  ISO<r  )  .GT.K  J>.  AND.  ITI  SOX  it)  .LT  .R  J1  )  )  GO  TO  4  C  G 
IF ( < TJSO<3) ,L T.F J).AND. FTISDI h) .GT .R J1 ) )  GO  TO  400 
IF  <  T I  SO ( b ) .EG  .RJ)  GO  TO  402 
GO  TO  snc 

400  CONTINUE 

>.  WRITC<7»401>  I.J.RAY(I.J.l) 

401  F0RMfcT(lX.*RAY(»»I3,I3»*.l)=»,Fb.2) 

CXLO  =  (K  AY<  I.  J»1  )-TISO(b>  >/<RAYt  l.J.D-R  A  YU,  <J*I)  .1  )  )*J 
EXCO=OS« (EXCD-l > 

COUNT  I 3)=COUNT(H)*I 
K=COUNT (B) 

EXRY<b»K)=EXCO 

EYrY(6»M)=DS*(I-I) 

CNTPT=CNTPT*l 
GO  TO  SCO 

402  CONTINUE 

COUNT <R)=COUNT<h)*I 
K  =  COUNT  «B) 
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EXf\  Y  <b*K  )-US*  (J-l) 

f  y  r  v  ( <i  t  k  >  -  r  :■  s  *  c  i  - 1  > 

C  N  T  P  T  =  C  N  T  p  T  *1 
b?C  CONTIMJE 

C  EXAMIfJt  TFMPS  VERTICALLY 
R  I  =  =t  AY  «  It  J,  H 
K I  1  =  K  A  Y ( < I -1 ) « J  t  1  ) 

IF  (  (HI  .GT.T  ISO{  1  >  >  .  AND.  <R  II  .GT.  T  ISO<  II  I  >  GO  TO  SF‘. 

IF(  (HI  .LT  .  T  IjC(MSO)  )  .  *\D.  (ft  I  1  .IT  .  TI  EO<  MSO)  >  )  00  TO  O'*. 

JO  b  2  b  n -  1  , N I  s ' ) 

IF  (  (  TI  GO  Cal  .GT.F  I  )  .AM}.  (TI  S0(  '• )  .LT  .R  1 1  >  >  GO  TO  bbZ 
I  F<  (  T  I  SJ<h  >  ,L  T  .  Kl  >  .  AMO.  (  T  1 .0  0  <  r  )  .GT  .K  II  )  )  GO  TO  tc>" 

00  TO  Mb 
bbO  COOT  I  V  J I 

C  ^IU(7,01)  I»J*-'AY(I«J*1) 

CYC  :  =  (FAY<<I-l)tJ*l>-TIS(T<‘.  >>/CFAY(tI-l>tJ*l>-K  A  Y(I,J«1)>»(I-1) 

i  rr.hz'j::  meycu-i  > 

CUlI'.TI  x  l=C  J'J\T  (  :  )  *1 
*=l GJ\T (H) 
tYiYlfiKUlYCO 
f.X,-Yli  .  r.  >  =00*  (  J-l  ) 

CNTPT  =  C>iTPTM 
b-’S  CO.nNJL 
CCO  COOT  I  \  J  f 
i :  o  co*  t i \ ut 

C  ACj  LOOP  IS  for  A1...HT  HA\j  SI.:; 

00  ?Co  I  =  Y  T  t  Y  « 

00  Gl  r-1 .  «i Of 
A  1  =  •>  A  Y  <  1  ,  X  •  1  > 

K I i -RA  Y  t  ( I -1 ) . X  ,  1  ) 

I  F  (  MI  .GT  .  <  TI  VJ  (0  )  )  >  .  AUu.t"  I  l  .LT.  (  T  I  SO  < >  )  >  )  GO  TO  fib's 
I  M  (R  I  .LT  .  (  T  1  V)(t< )  )  )  .  A\D.  (  x  1 1  .'T  .  (  TISC  l!  )  )  >  >  GO  TO  »b? 

IF  (  -  I  .Co,  .T  I  bO  (  >  )  G-.-  TO  F  b  * 

GO  T  0  p b 1 

a  b  z  c  o  M I  r«  j  r 

•: yc j  =  <  i -n ♦<•<  1 1 -t  is  j<i>) » /<.< n  -  ay<  i.  x«  n  > 

CYC.  =  G  *  <  F  Y  C  ">  - : » 

COuM  <  j  T-G  JU'x  T  <••  )  ♦  1 
k=c./J*tt  <  .) 

(.  YF  Y  <p,-<  )rf  YCo 
CXF  Y( Of  •  )  =  0.s»  (  x-1  > 

C  M  H  T  =  C  b  T  P  T  ♦  1 
GO  TO  £*3  1 
CGMIMJf 

CO- 'll  <!  ):nuU(  )»1 

<-C JO  Nil < !  » 

E  Y»  Y  (dtK  >  =  <  1  -1  > 

t  XA  Y  U't K  )  =  :■£»(  X-1  ) 

C  %  T  F  T  =  r  N  T  P  T  ♦  1 
tib  1  COOT  IMJl 
AGO  CO  »r  INJt 

C  M>';  LPoP  I G  T  '  H  TOP  GO* 

XX=X-1 

1=1 

JO  b  b  C  J  =  X  L  *  X  R 
R J  =  p AY  T 1 » J* 1> 

K  J1  =h  AY  <  1  »  <  1  >  ♦  n 

DO  bul  OsltVlSO 

I F C < TlGOCd > .GT.- J)  .-M). (T I '  0 C ■  >  .L T  .  Jll)  GO  TO  *<b 
IF  (  (T  ISoit-T  .LT  .*J>  .AND.  <  Tib  J<  l.,T,-  ot  »  >  '■ TO  •</,; 

I  F  1  F  j  .  F  .  T  l  i  0  (  )  )  G 1  -  T  O  h  b  ' 

GO  TO  Tbl 
bbC  COM  INJt 

tXC'J  =  <KAY(I,J,ll-Tl:0{:>)/<'-AYJI,J,l)-w,Y(J,tJ*)),  H  >  ♦  J 

cxr.  -zi)  , * <  c  xco-i  > 

COUNTY  -1  )  =CO  UG  T  (  - » ♦; 

K  =  C  jUN  T lu> 

CXPY(d,K)=CXC J 
EYI,  V»JtK>  =  '. 

CNTPT=CNTPT*1 
GO  TO  Hf.l 
H6i  COM  IMJl 

COOriT<P»=CGUM«  )«1 
R=Cvj'JMT  (i  ) 

EXRY<d«i‘.)=JG*<J-l  > 

CYKY(H.KJ=0 
CMPT  =  C0TPT«1 
3bl  COVTIMJC 
HbO  COGTINJC 

«-(I  TF  «  7, RC>  CMPT 
DO  Cbb  o=l»NISO 

C  =C OUNT »  0) 

WRITE  (It  iitL 
WRITETTtdl)  b.L 
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91  FORMAT ( IX,  * 1S0THH  »*,I2,*  HAS*,I4»*  POINTS*) 

IF(L.E3.0>  GO  TO  f ,6b 

WR I  TEC  7*90)  (EXRY(B»K),EYRY(b»K)»<=l»L> 

90  FORMAT  <1X»2F15. 3) 
h6b  CONTINUE 

WRITE (1*93)  CNT  FT 

93  F ORMAT ( IX, • TOTAL  NO.  POINTS  F0UN0=*,I5> 

92  F0RKAT(  IX, 15) 

CALL  CONTRLO,  'POINTS*,  7) 

RETURN 

END 

C  ****»♦»***.»»»*»**»*»«»*«.*»**»»*♦*»»******»***»»,*** 

C  SSDATA,  DATA  FOR  SSCONDJCT 
SUBROUTINE  SSDATA 
IMPLICIT  INTEGER(A,j,Q,W,X,Y,Z> 

C0MM0N/M1I /  A,X,Y,NIS0 

COM MON /MIR/  0S.DI,TMAX,TM1,\,FM4),THK(2>,H(2).T1S0(9)# 
s.  R  A  Y  ( I  Z  C  «  1  C  0 * 3 ) 

ALL  UNITS  ARE  IN  ME T C R S ,H0 JR S , CE L S] US ,  PJT 
IF  YOU  JSE  ENGLISH  UNITS, MAKE  CURE  THAT 
THEY  ARE  CONSISTENT. 

D  S  =  D IS  T  A  N  C  E  BETWEEN  NODES  (M) 

TSRF=SURrACE  TEMP  (OUTSIDE  GRID)  (INFINITE  OIST  AWAY) 

TriTR-uOTTOM  TEMP  (UNDER  GRID)  (INFINITE  CIST  AWAY) 

TR1T-TEMP  TO  RIGHT  OF  GRID  (INFINITE  HIST) 

TLFT-TEMP  TO  LEFT  OF  «,KIO  (INFINITE  GIST  AWAY) 

X =  NO.  Or  GRID  NODES  HORIZONTALLY 
Y -  NO.  OF  GRIO  NODES  VERT1CALLLY 
A-  .NO.  OF  DIFFERENT  MATERIALS  IN  THE  GPU 

DI  =  DISTANCE  HALF.AY  TO  INFINITY  (DI  IS  THE  MULTIPLE  Or  PS) 

TH*(L)  =  THERMAL  CONDUCTIVITY  (J/M*K) 

H(L)=C0NVECTI0N  COEFFICIENT 
K  I  (  1 )  =  HEAT  FLJX  FROM  TOP  OF  GRIO 
F I ( 2 ) -  HEAT  FLJX  FROM  LEFT  SIDE  OF  GKIU 
F 1 ( 3 ) -  HEAT  FLJX  FROM  30TT0M  OF  GRID 
F  I  (  4  )  =  HEAT  FLJX  FROM  RIGHT  SIDE 
RAY(I,J,l)=PREStNT  NODAL  TEMP 
RAY<I,J,r>-LOCAT2&.N  TYPE: 

=  1  INTERIOR  NODE  ,  FLUCTUATING  TEMP 
=2  NODE  ON  CONSTANT  TEMP  BOUNDARY 
=3  ON  CONSTANT  FLUX  BOUNDARY 
=  4  ON  CONVECTIVE  SUxFACE  E-GUNPArY 
=  3  INTERIOR  NO  J  E  ,  CONSTANT  TEMP 
=  fi  St  MI-INFINITE  f.CUNOANY 

=  7  NUDE  ADJACENT  TO  LEFT  SEMI-INF  riNDR  V 
=  :•  NOOt  ADJACENT  TO  DTTOM  Si  MI -INF  r<NDRY 

z  ■(  N DUE  ADJACENT  TO  RIGHT  SEMI-INF  9NDRY 

=  ld  NCOt  ADJACENT  TO  TO3  SlMI-INF  PNDRY 

-11  CONSTANT  CORNER  NDLl 

=  1?  CONST  FLUX  (ON  :J  0  T  H  t  SICES)  CORNER  NOPT 
=17  CONNECT 1VE (ON  rOlH  SUES)  CONNER  NOPE 
=  14  SEM-INF  (ON  BOTH  SIDES)  CORNER  NODE 
=13  SQUARE  NODE  ADJACENT  TO  TWO  SEMI-INF  SIPES 
=  3  C  SEMI-INF  NO  JE  ABOVE  SEMl-JNF  CORNER  NODE 

=31  SEMI-INF  NOuE  TO  LEFT  OF  SEMI-INF  CORNER  NODI 

=32  SCMI-INf  NODE  hL LOW  StMl-INF  CORNER  NODE 

=73  SEMI-INF  NODE  TO  RIGHT  OF  SEMI-INF  CORNER  NODI 

=  10  CORNER  » I T  H  VERTICAL  SITE  C ONV ECT , HOR I L  SIDE  CONST  FLUX 

=17  CORNER  dITH  HOAIZ  SIDE  CONVFCT.VTRT  SIDE  CONST  FLUX 
= 1 H  BOTTOM  SIGHT  CORNER.  VERT=CONST  FL UX ,HOR I Z = SEM I -1 NF 
=  lr-  MC-HT  SIDE  CONST  FLUX  NOPE  ADJACENT  TO  BOTTOM  STMI-INF 
=  2u  TOP  LEFT  CORNER.  Vf K T  =  SE Ml  -  1 NF , HOR I Z  =  CONST  FLUX 
=21  TOP  CONST  FLJX  N03E  ADJACENT  TO  LEFT  SEMI-INF 
=22  BOTTOM  LEFT  COENEr.  VERT=C3NST  FL JX .HOR I Z=SEMI -I NF 
=23  LEFT  SIDE  CONST  FLUX  ADJACENT  TO  BOTTOM  SEMI-INF 
=  24  UTTUN  rlGHT  CORNER.  V F R T  =  CO N VEC T ,  HOE  1 Z  =  SE“ I - I NF 
=  23  RIGHT  SUE  CONVCCT  NODE  ADJAC  TO  BOTTOM  SEMI-INF 
=  2b  TOP  LEFT  CGRNER.  V  EP.T  =  Sf  MI  - 1  NF  ,H  OR  I  Z  =  CON  VE  CT 
=27  TCP  CONNECT  NODE  ADJAC  TO  LfFT  SEMI-INF 
=  2-  i*  OT  T  DM  LEFT  CORNER.  VE  R  T  =CON  VEC  T  ,  HOR  I  Z  =  SEM  I  - 1  NF 
=  ?9  LEFT  SUC  CONVECT  NOUE  ADJACENT  TO  BOTTOM  SEMI-INF 
RAY(I,J,3)a  INDEX  OF  MATERIAL  (RANGES  FROM  1  TO  A) 

N I  S  0  =  NO.  OF  ISOTHERMS  TO  t)E  PLOTTED 
TMA<  =  TEMP  OF  HOTTEST  ISOTHERM  (*C) 

TM [ N=  TEMP  OF  COLDEST  ISCTHEhM  (»C) 

TISO(t)=  ARRAY  CONTAINING  ISOTHERM  TEMPS 
TISO(l)  HAS  THE  HOTTEST  ISOTHFRM 
COUNT <3)=C0UNTER  FOR  NO.ELTS  IN  EACH  ISOTHERM 

INITIALIZE  VARIABLES 
0S=  .  1 
A  =  1 
Xi?I 
Y  =  31 
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0 1  -  SO • 

TSRF  .  .  .  T:,TM  ONLY  MATTER  FOX  CONST  OX  CONVECTIVE  OK  SEMI -IMF  F.OL'NDARY 
F  Ok  CONNECTIVE  '.\IDK  Y  »  T  SRF  .  .  .  T  X  T  M  INDICATES  Tt.MX  A  V  A  Y  FROM  OKU). 

F  OK  COME  T  AMT  bOiJNOKY  «TSKF  .  .  .T  II  X  INdCATL  THE  TtMP  OF  Inf  hOUNPAR  Y . 

FOR  SEMI-INF  uuUNDHY  »TSRF  .  .  .1  f.TM  lNLlCATr  Tnf  TtMP  AT  INFINITY. 

tskfuc.o 

TLF  T-l 0  * 0 
TRITUG. 0 
TDTMUO.O 

THE  FJLLObllM.  VALULS  <ATTIX  '’ML  Y  IF  YOU  1MTINP  TO  LOCATE 
AN  J  PLOT  THL  ISOTHLkMS  OF  THE  RESULTING  TCMPERATJRr 
0  I  S  T  R  I  :j  U  T  1  0  M  .  IF  YOU  00  MOT  Vl'jH  TO  RUM  0(j"ROJTlMF 

IbOTHrt,  LEAVE  THE  FOLLOWING  ••  L  IMS  A'.  They  AFF. 

T  *U  K  z  1  3  0 . 0 
TMlNUn 
N  I  S  0  I A 

IM1  LALl/l  ISOTHERM  TEMPO 

l  Nl  ORDER.  wITH  TIEO(l)  HOTTEST 

m;.«i>  =  ioc.a 

TIS,-,(2l:7i/.0 
T  lcM.><i)rOO.O 
T I  SIX  A) -Si)  .  0 

IMTIALIZI  FLUXES  tKO“  u  0  U  M  i)  A  i .  I  c  r. : 

_,l  T  TO  ZERO  IF  MOT  USE  OF  IF  1  FOUL  '  T  E  0 , 

F  I  <  i  )  :  J  . 

f  i  < ;  >  =  j . 

Pl<_')-f. 

nio:;. 


!) j  MOT  CHAMGl  lor  fgllo-img  u  limes 
UJ  ■’  UU»X 
P  0  *j  I  —  1  *  Y" 

Li  J  /  X  1  t 

RAVI  [ t  J  *  K ) - ; 

7  CONTINUE 

b  COMTIMoP 

C  0 , 1  I  Nl  J  E 
J  J  =  I  •  A 
T  n  «,  <  >  )  :  " 

F'  <  )  - 

h  C  Ci.  T  i  M  u  - 


SFT  UP  •1(1, J,.')  1  t  •  OF  p  A  T  '  a  I  A  L  j 

T  t4 L  a a R  u  A  i .  A  T  »  a  I  A  L  .  •  L  ‘  E  To*  UJuN  Oi/’F  I  ‘j 

•'ATE  a  ML  l  E  H  F  *  /  l  -i  *  •■’ATI  a  I,.  L.  FT  AaT 

.Hi  U*. 

I  .  M  C  a  T  t  aAY<  I  *.!,')  r  A  ETC-  I  A  i. ''  ■<  1  TO  MHII 
this  le  ap  pssi.us  the  »mh*  *»amk1‘L  r  ; 

T  H  t  P  ;  A  l  \  I  N  >  IF  i  ii  .0  *  -i  J  -  H  A  0  w  A  T  F< ;  E  I X  • 

."O  1,:  JMO 
l.  J  A  I  r  1  »  Y 

IF  MAY  «I  ♦  J  t  *)  .\t  .;  >  1.-.'  TO  . 

X  A  Y  (  1  ,  J  ,  .5  )  I  A 
2  0  CV.  UN  Ji¬ 
lt)  C  0  v  T  I  M  j  l 

set  o°  material  pkcuktus  f  u  t  ai:<  ''atixial. 

THf  <  A) U . Aa 
M  <  A  >  -  0 


SET  JP  RAY  <1  i  J  *  A)  \ 0 j  x  L 

F’OjNOAR  YS 

CHANGE  ONLY  T  FT  f  "  R  Y  (  I  ,  J«  2 )  = 
Yf- Y-l 


C  TOP  hCJVOARY 
DO  AO  J  =  1,X 
RA  IIUJ.2)'-? 
KAYU*J«1)=TSKF 
AQ  CONTINUE 
C  E'OTTOM  t'OUNUAKY 
00  SO  j:|,< 

R A Y I Y» J.2) =u 
RAYTY* J«l)=TuT  M 
SO  CONTINUE 
C  LFFT  E’OUNPARY 
00  63  I=2»YT 
RAr«l,l»2)=6 
RAY(Itl»l)iTLFT 
60  CONTINUE 
C  RIGHT  UOJNOARY 


LOCATION  TYPE 
»  STM:  HIM. 


DO  7y  1=2. YT 

ra y  < l»  x.2»  =.' 

72  CONTINUE 

CORNERS 

AJJjST  fOTh  RAY(I,U.2J  4  Ml  RAYtl.U.ll 
7  OP  LEFT  CORSE R 
RAY  (1.1.0  =  11 
RAY(i.i.n  =  io.c: 

E0T1CM  LtFT  CORNER 
R AY ( Y • 1 .2)  =  14 
R AY  x Y, 1. 1)  =  10 .0 
lOTTOM  rIGhT  CJrMER 
RAYJY«X«2)  =  l.-> 

R4Y<Y,X,1»=1C.C 
T 0°  RIGHT  CORNER 
s  4Y< 1, x.21  =  11 

RAYd.x,  11  =  10.  r 

INTEaIOR  NCdES  ?  I  Vl  \  FLuCIUAI  IMG  UMP 
2  0  N  •  T  CHANGE  THIS  L  00°  • 

Y Y=Y -1 

>  x  =  x  - 1 

c  o  ;  j = o .  x  x 
o  o  - :  i  =  2 .  v  Y 
RAY<I.J«0)=; 

’  continue 

:  couiw 

IU'ICATl  LOCATION  of  all  imomok  nodes  of  const  TYMP 
Cl  T  RAY<I.J»2)  =  C  FOR  LOCATION  T  Y  F  £  1 
-v  AT  <  10,  X.2  )  =  E 
R  4  Y  X  1 1  ,  A  ,  2  )  =C 
xxriiii ( x - l >,2)=s 
•\Av<12,X»2)=d 

T *t  FOLlO.ING  IS  FOR  Lf  r  T  SIDE  X  BOTTOM  SI  MI -IMF  I  MITT. 
A2JCST  T-ESE  LOCFS  AM'  NUMBERS  TO  FROVl'f  FOR 
SOLES  ADJACENT  TO  St M I -I N1 I M I TE  3OJN0ARY  MOOES. 

YT=Y-l 

OC  ->1  I  =  2  «  Y  Y 
R  A  T  (  i  *  2 . 2 »  =  7 

1  CONTINUE 
x  x  =  X -1 

20  -2  J  =  2  »  »  X 
RAY(«Y-3),J,2>=- 

2  CONTIMJF 

R A  r ( I Y -1 >  *  1 .2  >  =0  1 
RAY( <Y-1),X,2)=19 
RAY»(Y-l),2«i)=io 

SET  OP  -AY(I.J.l)  ...MODAL  TEMPERATURES 
IF  CONuT  T"«P  rOJSDARYS  ARE  SOT  UNIFORM  TC'P, 

(TSxF.T  TM.TLFT  .TrvIT  ,  «  r  <'-F  E  C  T  IVELY  )  .THEN  INDICATE 
TEMPERATURES  OF  sOJMDARY  MOOES  HERE.  ALSO  INDICATE 
TEMPERATURES  OF  CONSTANT  INTERIOR  MODES  HERE. 
RAYm,x.l>=lC0.0 
R  AY  ill » A,1 ) =1C0.3 
RAY  <  H , ( X- 1 ) ,  1) =1CC  .C 
RAY(12,Xfl)=130.0 

JO  MOT  Ch-MOl  THE  FOLLOWING  If  LINES. 
•RlTEtD.lil)  JS.DI 
SI  PDM',AT<lX,2F'>.<t) 

»RITE<3,1J2)  A.x.Y.NISU 
12  F  G  s  “  A  T (1X.4I3I 

.Slut  j,lJi)(TRML).L--l»AJ 
•xIT2<C,133>  <n(LJ,L=l,A» 

.R1TE15.13J)  <TISC(3>»:R=1,MS0I 
Jo  FOSMAT I1X»F9«<»J 

■nlTECS.liA)  (FI(I),I=1,4> 

3<*  FOR'ATC 1X..F9.A) 

,R1TE(2,151>  (<RAY(I»J,3),J=1,X),I=1»Y> 

«RITl<5»15)1>  <<rAY<I,U»2)»J=1.X)»I=1*Y) 
yRITElb.lSl)  (CAAY<I,J,1),J=1,X>,I=1,Y) 
bl  F0R“AT«1X,11F7.2» 

RET  JRM 
ENJ 


d  A  S  0  m  X 
USAGE 

ARGUMENTS  A 
N 


CALL  8  AND  MX  ( A  »M  »NLC  »NUC , I A ,8 ,M , 18, I JOB. XL. 
IERJ 

1 NPUT /OUTPUT  MATRIX  OF  DIMENSION  N  BY 
<  MUC  *.ML  C*1 )  .  SEE  PARAMETER  IJOB. 

ORDER  OF  MATRIX  A  AND  THE  NUMBER  OF  ROWS  IN 
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NLt  "  NUMBER  ^  GF^LOWEK  COD  1  AGON  AL  S  IN  MATRIX  A. 

\UC  -  N  U  t  f  ft  ^  0  K  UPPER  COD  J  Ao  ON ALS  IN  MATRIX  A. 

11  _  ft  o»  *  DI  ”C  V.  1  ON  OF  MATRIX  A  EXACTLY  AS 

''  SPEC 1M l U  IN  THl  DIMENSION  STATEMENT  IN  THE 

CALLING  PROGRAM.  (INPUT)  „  _ 

u  -  IN  PUT  / OUTPUT  MATRIX  OF  DIMENSION  N  *  Y  M. 

ON  INPUT*  0  CONTAINS  THE  M  R1 CHI -HAND  S IDES 
OF  THE  EQUATION  AX  =  ".  ON  OUTPUT.  THT_ 
SOLUTION  MATRIX  X  REPLACES  r-..  IF  1JO.  -  1. 

M  -  NUMBER  OF^IGHT ‘HAND  SIDES  (COLUMNS  IN  H>. 

in  -  Koi Dimension  of  matrix  b  exactly  as 

SPECIFIED  in  THE  JIMEVSION  STATEMENT  IN  THE 

-  -  «t  iKHW' 

EQUATION  AX  =  B .  ON  INPUT.  A  CONTAINS  THE 
COEFFICIENT  MATRIX  OF  THE  EQUATION  AX  =  B. 
UHCRF  A  IS  ASSUMED  TO  Bf  AN  N  BY  N  BAND 
MATRIX.  A  IS  STORED  IN  BAND  ST0RA6C  MCOE 
AND  THEREFORE  HAS  DIMENSION  N  BY 
( N  L  C  ♦  iNU  C  ♦  1  )  »  ON  OUTPUT.  A  IS  REPLACED 
BY  THE  U  MATRIX  OF  THE  L-U  OECOMPO S IT  I  ON 
OF  A  ROWWISE  PERMUTATION  OF  MATRIX  A.  U 
» c  STORED  IN  3 ANU  STORAGE  "ODE. 

i  =  •,  factor  the  matrix  a.  a  contains  THt 

SAME  INPUT/OUTPUT  INFORMATION  AS  IF 

T  SOLVE  THE  EQUATION  AX  =  B.  THIS 

OPTION  IMPLIES  THAT  BANDMX  ALREADY 

pprn  CALLED  USING  IJOB  :  !  OR  1  SO  THAT 
THE  MATRIX  AHAS  ALREADY  BEEN  FACTORED. 

IN  THIS  CASE*  OUTPUT  MATRICES  A  AND  XL 
MUST  HAVE  BEEN  SAVED  FOR  REUSE  IN  THE 

VI  -  u  OHKCAft  E  A^  OE^D  I  MENS  I  ON  N*(NLC-»tl.  THE  FIRST 

*L  NLC*N  LOCATIONS 3F  XL  CONTAIN  COMPONENTS  OF 

THE  L  MATRIX  OF  THE  ^DECOMPOSITION  OF  A 
ROWWISE  PERMUTATION  OF  A.  THE  LAST  N 
LOCATIONS  CONTAIN  THT  PIVOT  INDICES. 

IEK  -  ERROR  PARAMETER.  (OUTPUT) 

TE I ERN=L ll^INDI CATES  THAT  MATRIX  A  IS 
L algor  I THMI CALL Y  SINGULAR.  (SEE  THE 
CHAPTER  L  PRELUDE). 

REQUIRED  SUbROUTINFS-SUBl.SU _ 

SUBROUTINE  BANDMX  ( N.NLC .NUC ♦ I A »M. I E . I  JOB.  1  ER ) 

COMMON/ MO/  A(lOCOA,?Ob>.P«lPu0A»b).XLCUCUA.lO5) 

DIMENSION  A ( Y » X.NUC *NLC* 1 )  ?“D 

XL ( Y«X  .  (NLC*1 > >  1-U 

B(Y.X.l)  P- D 

DArA  ^EKO/0'/,OnRST°EXECUTABLE  STATEMENT 


IER  =  0 
JBEG  =  NLC  ♦  1 
NLC1  =  UBEG 
IF  (IJOB  .EQ. 
RN  =  N 


2)  GO  TO  FO 


RE  STRUCTURE  THE  MATRIX 

FIND  RECIPROCAL  OF  THE  LARGEST 

AtSOLUTE  VALUE  IN  ROW  I 


I  =  1 

NC  =  JBE  G*NUC 
NN  =  NC 

fp^IN  .EU.  1  .OR.  NLC  »EQ.  0)  GO  TO  C5 
5  K  =  1 
P  =  ZERO 

DO  10  J  =  JUEG.JEND 

A  ( I  »K  )  =  A ( I »U> 

Q  =  ABS ( A ( I  *  K  >  > 

IF  (0  . GT.  P)  P  =  w 
K  =  K  ♦! 

10  ii.no>  00  to  135 

Jl'Il’iS,”  5cf"sOPtO  20 
00  15,M  5*3U 


i 


15  CONTINUE 
20  I  =  1*1 

I ( JEND-JBEG  .£Q.  N>  JEND  =  UEND-1 
IF  (I  .LC.  NLC)  GO  TO  5 
UUE  3  =  I 
NN  =  JEND 
25  JE NO  =  N-NUC 

DO  1  :  J6EG*N 
P  =  ZERO 
DO  30  J  =  1 »  NN 

Q  =  ABS(A(I*J>) 

IF  <a  .GT.  P)  P  =  U 
f  ON  T INUE 

-  ZERO)  GO  TO  135 
=  ONE/P 

JENO)  GO  TO  37 
JENO) 


3  C 


35 

37 


if  (P  .Ea. 

XL ( I *NLC  1 ) 

IF  <1  .Ea. 

IF  <1  .LT. 

<  =  NN*1 
00  35  J  =  K ♦ NC 
A  (  1  *  J )  =  ZERO 
CONTINUE 
NN  =  NN  *  1 


GO  TO  AO 


AO  CONTINUE 


L 

00 


NLC 


L-U  DECOMPOSITION 


45 

50 


75  K  =  1 »N 

P  =  ASS < A <Kt 1 ) ) »XL <K «NLC1  ) 

I  =  K 

IF  ( L  .LT.  N>  L  =  L ♦ 1 
<1  =  K ♦ 1 

IF  ( K 1  . GT •  L)  GO  TO  51 
00  a5  J  =  K 1  »  L 

a  =  ABS  (  A ( J« 1> > *XL ( J »NLC1 > 
IF  <Q  .LE.  P>  GO  TO  A5 
P  =  3 
I  =  J 
CONTINUE 

XLUtMCl)  =  XHK.NLCl) 
XL«K,NLC1)  -  1 


=  RN*P 


I- 

(G 

.EQ 

IF 

<K 

•  E  U 

r  0 

55 

J  - 

RN )  GO  TO  1 3i 

,  1)  GO  TO  60 
1«»'C 


SINGULARITY  FOUND 
INTERCHANGE  POJS  1  AND  K 


TO  75 


A  <K  » J )  =  A(I«J) 

A< I«J>  =  P 
CONTINUE 

IF  ( K 1  .GT.  L»  GO 
JO  7 J  I  =  K1  .L 

P  =  A(I«1)/A(K«X) 

IK  =  l-K 
XL(Kl.IK)  =  P 

JO  65  J  =  2  *  N  C  ^  .. 

A<I»J-1)  :  A  (,I  «  J) -P*  A  <  K  t  J  ) 
CONTINUE 
A  < I  *  NC  I  =  ZERO 
70  CONTINUE 

»'  ■>  GO  TO  000 


55 

Si 


63 


FORWARD  SUBSTITUTION 


BO 


95 

90 


99 

IOC 

105 


L  =  NLC 

DO  10  5  K  i  1  *  N 

I  =  XL<K«NLC1I 
IF  < I  .EQ.  X  >  GO  TO 
30  85  J  : 

P  =  0<K,J) 
h ( K  « J )  =  »i  I »  U ) 
t$  1 1 1 J )  =  P 
CONTINUE 

IF  IL  .LT.  NT  L  - 
Al  =  K ♦  1 

IF  (K 1  .GT .  L  >  GO 
DO  100  I  =  K1»L 
IK  :  I-K 
P  =  XL(K1»IKI 
DO  95  J  =  1 »M 

DlIiJ)  =  H  « 1  ,J)-P»B«K«J> 

CONTINUE 

CONTINUE 

C0N1INUE  h ACK  W  AR  0  SUBSTITUTION 


90 


L*  1 

TO  10‘ 
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<  innnnoononont  uc  ;o»  JOfKKioonnnono 


JBE  G  =  NUONLC 
00  1 2b  U  =  1 »  M 
L  -  1 
K 1  =  N*1 
00  120  I  =  1  ,  N 
K  =  Kl-I 
P  =  3<K,J> 

IF  < L  .E0.  1)  GO  TO  lib 
DO  110  KK  =  2  »  L 

j  u  —  ft  u  m 

p  =’p-A (K.KK) *b< IK-1  ,J> 
110  CONTINUE 

lib  B  {  K  »  J  )  =  P/ A  <  <» 1 ) 

IF  <L  .LE.  JBFG)  L  =  L*1 
120  CONTINJE 
12b  CONTINUE 

GO  TO  000b 
lib  I  EH  =  121 
1000  CONTINUE 

a  R  I  T  E  (  6  «  998  > 

198  FORMAT ( IX, ‘CALL  SUP 1 • ) 

CALL  SUE  1  < I FR  » • LE  *T  IB  •) 

900b  KlTURN 
END 


SUP1 


PJRP0SE  -  POINT  A  MES  SAGE  REFLECTING  AN  ERROR  CONDITION 

USAGE  -  CALL  SU91  I1ER. NAVI) 

ARGUMENTS  If-  -  ERROR  PAHAMLTLH.  (INPUT) 

ILh  :  1*J  WHERE 

i  =  126  implies  terminal  err^r, 

I  =  feA  IMPLIES  WARNING  «  I  T  H  FIX,  AND 
i  =  32  IMPLIES  WARNING. 

J  =  ERROR  CD JE  RELEVANT  TO  CALLING 
ROUI  1NE . 

NAME  -  A  SIX  CHARACTER  LITERAL  STRING  GIVING  THE 
NAME  GF  THE  CALLIN'-  ROUTINE.  (INPUT) 

REuUIRlD  SJi'RTlNF=  SUoL 

REMARKS  the  ERROR  MESSAGE  PRODUCED  UY  ^Ut>l  IS  WRITTEN 

ONTO  THE  STANDARD  OUTPJT  UNIT.  THE  OUTPUT  UNIT 
\  J  M '  E  R  can  .!,*  OElEKHI'lEr  -Y  CALLING  SU-?  t  3 
c0LL0«S..  CALL  SUED ( 1  .N1N.N0UT ) . 

T M  E  OUTPUT  UNIT  J  M  u  f  R  CAN  BE  CHANGE  0  PY  CALLING 
EUUL  AS  FOLLOWS.. 

m»  -  o 

NOUT  -  NEW  OUTPUT  UNIT  NUMBER 
CALL  3UL'M3,MN»N">UT  > 

SEE  THE  UGETlO  DOCUMENT  FOR  MOPE  DETAILS. 


SU'ROJTl.NE  SUB1  (IER.NAME) 

SPECIFICATIONS  FUR  'RGJWENTS 

INT'GER  1EH 

REAL*;)  NAME 

SPECIFICATIONS  FOR  LOCAL  VARIABLE 

kEAL*H  namsei.nameo 

OATA  NAMSF T/HHUrRSE T  / 

DATA  NAMEU/pH  / 

FIRST  EXECUTABLE  STATEMENT 
DATA  LEVEL /A /.I CuOF / u/ ,  I C  I/1H=/ 

IF  (IER.GT.909)  GO  TO  25 

IF  (IER.LT. -32)  GO  TO  bb 

IF  (IER.LE.12d)  GO  TO  E 

IF  (LEVEL. LT.l)  GO  TO  30 

PRINT  TERMINAL  MENS  A  G i 

CALL  SU«2( 1 »NIN* IQUN1T ) 

IF  CIE3UF.EQ.1)  UR  I  TE  (  I  OUNI  T  ,  3b  )  1  f  R  »N  AMI  Q  .  I  E ..  « »'A  *C 
IF  (IlJDF.EO.O)  UK  I TE  < I OUN I T  »  3b  )  IFR.NAIE 
GO  TO  3C 

IF  (IER.LE.faA)  GO  TO  1C 
IF  (LEVEL.LT. 2)  CO  TO  30 

print  warning  *.ith  fix  message 

CALL  SUB2T 1 .N IN, I OUNlT  ) 

IF  <ICv.OF.Eu.  1  >  UK  I  TE  (  IQUNI  T  ,  At  )  I  F  R  ,  \  t  *lt  »  ,  I  f NO  ,  .V  A  Mf 
IF  (IEUOF.EG.O)  W«I  TE  UGUM  T  ,  AD  >  IER»NAWF 
GO  TO  3C 

IF  (ICR. LE, 32)  GO  TO  lb 

PRINT  uARMNG  HESS  A  GF 

IF  (LEVEL.LT. 3)  GO  TO  30 
CALL  SUB2 ( 1 ,N  I  .N  ♦  I OUN IT) 

IF  (IEuOF.EU.l)  UR  I  TE  (I  OUM  T,  As)  U  R  ,N  A  MF.)  ,  1 1  U  ,  N  AMI 
IF  (IEuDF.EG.0)  U  R  I  T  E  (  I  OUN  IT  ,  A  b )  Ii'F.NAvf 


57 


ci  clary  a  r>oc;oc >t  n  jooc  mcu  ii  ioooi  jooooooc.ooo  r*on 


15 


25 


GO  TO  30 
CONTINUE 

IF  (NAME.NE.NAMSET)  GO 
LEV3LD  =  LEVEL 
LEVEL  =  IER 
IER  =  LEVOLO 
IF  (LEVEL. LT.O)  LEVEL  = 
IF  (LEVEL. GT. 4)  LEVEL  = 
GO  TO  30 
CONTINUE 
IF  (LEVCL.LT. 4) 


CHECK 
TO  25 


FOR  UERSf T  CALL 


GO  TO  3 C 

PRINT  NON-UEFINED  MESSAGE 

CALL  SUd2(l,NlN,IGUNIT> 

IF  (IEJOF.tO.l)  WHI  TE ( IOUNIT ,50 >  1 ER »N AMI Q *  1 EQ , NAME 
IF  (IEJDF.EO.O)  WR I T E (  I  CUN  IT  »  51  >  1  Ek  ,NAME 


TERMINAL  ERROR. 13X,7H(IER  =  .13. 
ROUTINE  .AG.Al.AG) 

kARNING  WITH  FIX  ERROR  (IER  =  ,13* 
ROUTINE  .A6.A1.AG) 
kARNING  ERROR, 11X,7H(IER  =  ,13, 
ROUTINE  .A6.A1.A6) 

UNDEFINED  ERROR, 9X,7H(IlK  =  ,15, 
ROUTINE  .A6.Al.Ab) 

SAVE  P  FOR  P  :  R  CASE 
P  IS  THE  PAGE  NAME 
R  IS  THE  ROUTINE  NA«E 


30 

IEQDF  =  0 
GO  TO  65 

35 

FORMAT(l°H 

**  * 

S20H)  FROM 

IMSL 

40 

FORMAT (36H 

*  *  * 

S.2CH)  FROM 

I  MS L. 

45 

FORMAT! I8H 

*  *  * 

&20h)  FROM 

IMSL 

3  0 

FORMAK20H 

*** 

L20H)  FROM 

IMSL 

55 

I  E  5  OF  =  1 

NAMEO  =  NAME 

65 

CONTINUE 

RETURN 

END 

SUH2 


PURPOSE.  -  TO  RETRIEVE  CURRENT  VALUES  AND  TO  SET  NEW 

VALUES  FOR  INPUT  A  NO  OUTPUT  UNIT 
IDENTIFIERS. 


USAGE 


CALL  SUH2 ( IOPT .NIN, NOUT > 


ARGUMENTS 


REMARKS 


IOPT  -  OPTION  PARAMETER.  (INPUT) 

IF  I  OPT  =  1  ,  THE  CURRENT  INPUT  AND  OUTPUT 
UNIT  IDENTIFIER  VALUES  ARF  RETURNED  IN  NIN 
AND  NO  J  T  ,  RESPECTIVELY. 

IE  I  C  P  T  -  (31  THE  INTERVAL  VALUE  OF 

NIN  (NOUT)  IS  RESET  FOR  SUBSEQUENT  USE. 

NIN  -  INPUT  UNIT  IDENTIFIER. 

OUTPUT  IF  I  OPT  =1 «  INPUT  IF  I0PT=2. 

NOUT  -  OUTPUT  UNIT  IDENTIFIER. 

OUTPUT  IF  I0PT=1,  INPUT  IF  I0PT=3. 

E  A  Cm  IMSL  ROUTINE  THAT  PP  RF  OK  MS  INPUT  AND/OR  OUTPUT 
OPERATIONS  CALLS  SUU2  TO  OBTAIN  TXE  CURRENT  UNIT 
IDENTIFIER  VALUES.  IF  SUB2  IS  CALLED  WITH  IOPT=D  OR  3 
NEW  UNIT  IDENTIFIER  VALUES  ARE  ESTABLISHED.  SUBSEQUENT 
INPUT/OUTPUT  IS  PERFORMED  ON  THE  NEW  UNITS. 


5 

1  C 
9005 


SUBROUTINE  5062 ( IOPT , NIN, NOUT  ) 

SPECIFICATIONS  FOR  ARGUMENTS 

IMPLICIT  1  NT  E  GEK  *  4 ( I -M ) 

IMPLICIT  DOUBLE  PRECISION(A-H.O-Z) 

INTEGER  IOPT, NIN, NOUT 

SPECIFICATIONS  FOR  LOCAL  VARIABLES 

INTEGER  NIN U, NO U TO 

DATA  NIN 0/5/ »N0 U 10/6/ 

FIFST  EXECUTABLE  STATEMENT 


IF  (IOPT.EW.5)  GO  TD  1! 

IF  (IOPT. EO. 2)  GO  TO  5 

IF  (IOPT.NE.l)  GO  T'U  9CC5 

NIN  =  NINO 

NOUT  =  NOUT  D 

GO  TO  9305 

NIN'.)  ;  NIN 

GO  TO  9305 

NOUTD  =  NOUT 

RETURN 

END 
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Table  A1.  Sample  output  from  SSCONDUCT  —  provides  initial  information  and  final  temperature  distribution. 
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REST  OF  THE  COLUMNS 
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IE*=12^t  M  AT  h  I  X  I  S  SI-.o'JL  • 

IC^=0*  everything  is  ORAY. 
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-4tOlT^ChO<MtOin  yflNCOO'ONOOHHHHH^HOOOOff'ff'O 
CL  — #-4^4  »-t^<  CM  CM  CM  CM  CM  CM  CM  CM  CM  tO  tO  tO  fO  tO  tO  tO  tO  tO  fO  fO  tO  tO  CM  CM  r4 
*. 

Uf 

H-  O  *■  X  CM  \D  tr»  lO  >D  0s  CM  +  ^  OOO  *-4  CM  CM  CM  CM  *H  CT  ®  IP  tO  O  >0  O'  ♦  O'  ♦ 
O  O  0 1-4 1-4  «-4  CM  CM  CM  tO  fO  fO  tO  4T  <T  ♦  ^  +  4T  tO  tO  tO  tO  tO  CM  CM  ^4  #-4  o  O 

J  ••*•••*•  . . .  C  •  I  •••••••  • 

4000000000000000  OOOOOOOOOOOOOOOO 


JO  J-OO'N  ^f^iO^lDCJ'^'H^fOCCih-  n<M  4^^/»sfl  t)H»04000 

oo ocm in *0^0^*  sDoNf-tooovO^HOH^nHj'XcooOiin^N 
^^•••••••••••••••■•■••••••••••••* 

OsACM  XinCM^tOCM  £OC0*O<r  0sinCMOXX*-t0*-4O0'C^  x^  x  Q 
uj  *-i  cm  cm  •o  *  in  vdc^  (no  xr~p-Mvnu‘>inm*’*  *  *■  *  *totototo«o-i 
X  f~* 

t— 

OMOaOOO^N'ff'tO  *  (MC'J**NtOoSvfl(V(J'<f  *®C\jNo(MHCO 
LO  WCOff'C^tOS  Xin^4h» CJ'sHtOK  an>£fT'iniAcDNff'h'43SaOH4cOCNj 

Oinr-4f^*-^4XvDin*-£T'X»-iXr4h-*-r4(T'h'in*-CM»40(TNX  X^  \0O 

h-f-4^4CMCMto*-  xxnnm  *■*•*•*•*■  <r  *»otO'0,o^o^H 

CO 
Ui 

aco^4«h*knr-«r  oooc  «-r  r^tOr4tof^*CMviJin<r  cm*- — «xx*o CT'CMCMf^ 
o  <x  cm  o  o  <m  ^  to  on  in  cmo  •-«  m  r*»  o  r^*  sa  ®  cm  ao  m  4  to  *  \A  «  cm  \0  cm 
^••••••••••••••••••••••••••••••« 

roi/iHf«*o(M/>cvco4NNinHcoin'Oocn\nn>orgHOfTNOo^  a>o 

r4  cm  cm  to  to  4  m  in  si)  io  vO  >a  tn  m  in  in  *>  *  *-  *  <r  *■  **  to  to  to  to  to  *4 
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Table  A2.  Sample  output  from  SSCONDUCT 
isotherm  locations. 


72 

ISOTHM  a 


ISOTHM  » 


ISOTHM  W 


ISOTHM  3 


HAS 

4 

POINTS 

1.900 

1.000 

2.000 

0.900 

2.000 

1.000 

2.000 

1.100 

HAS 

12 

POINTS 

1 .864 

0.800 

1.900 

0.779 

1.757 

0.900 

1.600 

0.852 

1.71ft 

1.000 

1.726 

1.100 

1.775 

1.200 

1  .R98 

1.300 

1.600 

1.231 

1.900 

1.301 

2.000 

0.752 

2.000 

1.326 

HAS 

27 

POINTS 

1.854 

G.6C0 

1.900 

0.587 

1  .645 

0.700 

1.700 

0.664 

1  .  e  0  0 

0.615 

1.557 

0.600 

1.600 

0.739 

1 .468 

0.900 

1 .500 

0.R55 

1.424 

1.000 

1.401 

1.100 

1.543 

1.200 

1.400 

1.119 

1.402 

1.300 

1.400 

1.276 

1.427 

1.400 

1.472 

1 .500 

1.642 

1 .600 

1  .600 

1.546 

1 .651 

1.700 

1 .600 

1.662 

1.700 

1  .  735 

1.674 

1 . 6  0  C 

1  .  ft  0  0 

1.781 

1.400 

1.1-07 

2.000 

0.577 

2.000 

1.916 

>  HAS 

j  ) 

POINT; 

1.515 

0.40C 

1  .  6  Q  C 

0,  "><t 

1  .700 

0  .54  i 

1  .  *  j  0 

0.35. 

1  .  i  0  0 

0.321 

1.076 

0.500 

1.500 

0.489 

1.400 

0  •  4  A  4 

1 . 5  (■  C 

0.405 

1.101 

0  .  .6  0  0 

1.200 

0.542 

0 . 9  0  5 

0.  700 

1.000 

0  .  6  9 

1.100 

0.6  01 

0.824 

O.eCO 

3  .  0  0 

0.743 

0  .  7  J  li 

0.700 

0  .  -  00 

C.H22 

0.594 

1.000 

0.630 

0.994 

C.700 

0.4C6 

0.484 

1.100 

0.500 

1.087 

0.3  75 

1.200 

0.400 

1.160 

0.264 

1.300 

0.300 

1.271 

0.200 

1.361 

2.000 

0.317 
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Table  A3.  Sample  input  to  SSCONDUCT  —  file  generated  by  SSDATA  (subroutine  of  SSCONDUCT). 


0.1000  50.0000 

1  21  31  4 

1.4400 
0.0000 
100.0003 
70.0000 
50.0000 
30.000C 


0.0300 

3.0  0  03  0 

.0  30  0 

0.0000 

l.CC 

1.00 

1.00 

1.00 

1.C0 

1 

.00 

1 .30 

1.0C 

1.0  0 

1.3  0 

1.00 

1.00 

1 

.00 

1  .  C  3 

1.03 

1 .00 

1.0  0 

1.0  0 

1 . 0  c 

1.00 

l.CO 

1.30 

1.00 

1 

.00 

1.33 

1.00 

1.0  0 

1.00 

1.00 

1.00 

1.03 

1.00 

1.00 

1.03 

1 

.00 

1  .33 

1.0C 

1.03 

1.00 

1.00 

1  .  C  0 

1.03 

l.CO 

1.00 

1.00 

1 

.00 

1.30 

1.03 

1.03 

l.CC 

1.33 

l.C? 

l.CC 

1.00 

1.00 

1.00 

1 

.  0  j 

1.03 

1.03 

1.0C 

1.00 

1.3  3 

i.c : 

l.ro 

1  .  u  0 

l.OP 

1.00 

1 

.33 

1.33 

1 . 0  c 

1.03 

1.30 

1.00 

1.00 

1.00 

1.00 

l.CO 

1 .00 

1 

.30 

1  .0. 

1 .00 

1.0  3 

l.CC 

1.3  3 

i .  c  r 

1  .00 

l.ro 

1.33 

1.00 

1 

.00 

1  .03 

1.00 

1  .03 

l . ;  r, 

1.03 

1.00 

1.03 

l.ro 

l.PO 

1.00 

1 

.03 

1.03 

1.0  3 

i .  n  c 

1.33 

1.3  3 

i .  3 : 

1  .  i.  0 

l.CC 

In  C 

1.00 

1 

.00 

1.03 

l.CC 

1 . 0  c 

1  .C  3 

i.:: 

1  .DC 

1.00 

l.CO 

1.00 

1.00 

1 

.30 

1.33 

1  .  Ci  ci 

*  .  u  3 

1  .-  C 

1.3' 

i.:c 

i .  r,  o 

l.CO 

1 .00 

1.00 

1 

,00 

1.03 

1 . 0  3 

1.03 

i : 

l.CC 

l.CC 

i .  r  1. 

l.CO 

1.00 

1.00 

1 

.00 

1  •  0  i 

1 . 3  c 

1.3^- 

l.CC 

1  •  1  J 

l .  r  3 

1  .CO 

l.CO 

1.30 

1.00 

1 

.00 

1  .33 

i  .or 

1  .  U  J 

1  .uO 

1.00 

1.00 

1.00 

l.PO 

1.00 

1.03 

1 

.00 

1.33 

1.0 

1.0  0 

l.CC 

1.3  0 

1  .  r.  3 

1  •  hi  0 

1.00 

1.30 

1.00 

1 

.00 

1.33 

1.03 

1.0  0 

l.CC 

1.03 

1  .  *»  o 

1  .00 

1.00 

1  .  P  c 

1.00 

1 

.00 

1.03 

l.CC 

1  .  U 

1  .  u  C 

1.33 

l.CC 

1.33 

l.ro 

l.CC 

1.00 

1 

.00 

1.3  3 

1.3  3 

1.00 

1  .33 

1. 30 

1.30 

1.00 

1.0  0 

1.03 

1.00 

1 

.30 

1.3  3 

1.0c. 

1.00 

1  .  '•  C 

1.3  3 

l .  o  r 

1  .03 

1.00 

l.OP 

1.33 

1 

.00 

1  .  t  'J 

i  i  0  o 

1  .Co 

1  ..  3 

1.33 

1.00 

1.00 

1.0  8 

1.03 

1.03 

1 

.33 

1.33 

i  •  o 

l.CC 

l.CC 

1.3  3 

1.00 

1  .03 

l.r.o 

1.00 

1.3  3 

1 

.00 

1  .03 

l.Ot 

1.33 

l.iiO 

1.  JO 

1.00 

1.00 

l.CO 

1.00 

1.00 

1 

.03 

1.03 

1  .  G  j 

1.0c 

1  .  C  Ci 

1.33 

1.00 

1.00 

1.00 

1.00 

1.03 

1 

.30 

1.33 

1.30 

1.03 

1.33 

l.C  3 

l.vC 

1.C3 

1.0  3 

1.03 

1.03 

1 

.00 

1.33 

1.0  0 

l.C  0 

1.3  0 

1  .  3  3 

1.00 

1  .  00 

l.ro 

l.OP 

1.00 

1 

.00 

1.03 

X  .03 

1.00 

1.00 

1.00 

l.OP 

1.00 

1.00 

1.00 

1.03 

1 

•  33 

1 . 0  h! 

1.0  0 

1.00 

1.3  3 

1.00 

1.00 

1  .PO 

l.PO 

1.00 

1.03 

1 

.33 

1.30 

l  .  C  u 

1.00 

1.00 

1.00 

1.33 

l.ro 

l.CO 

1.00 

1.0C 

1 

.03 

1  .  Ci  J 

1.33 

1.03 

1  .  C  0 

1.3  3 

l.CC 

1.00 

l.ro 

l.CC 

1.03 

1 

.00 

1  .00 

1.03 

I. 00 

1.0  3 

1.  JO 

l.OP 

1.03 

1.00 

1  .  J  0 

1.00 

1 

.00 

1.3-1 

i.0  0 

1.00 

1.30 

1.03 

1.0  0 

i  .no 

l.CO 

1.00 

1.03 

1 

.33 

1.33 

1.00 

1.00 

1  .C  3 

1.00 

l.C  9 

l.CO 

l.CO 

1.00 

1.00 

1 

.00 

1.33 

1.33 

1.00 

l.CC 

1.30 

1.0G 

1.00 

l.CO 

1.00 

1.00 

1 

.03 

1.03 

1.03 

1.0  0 

1.00 

1.00 

1.00 

1.00 

1.00 

l.CC 

1.30 

1 

.00 

1.33 

1.0  0 

1.00 

1  .  V  0 

1.00 

1  .  C  3 

l.CO 

1.00 

1.00 

1.03 

1 

.03 

1.00 

1.30 

1.00 

1.30 

l.C  3 

1.00 

1.00 

l.CO 

1.00 

1.00 

1 

.00 

1.03 

1.00 

1.00 

1  .cG 

1.03 

1.0  (i 

1.00 

l.ro 

l  .on 

1.33 

1 

.00 

1.03 

1.0C 

1.00 

1.03 

1.00 

1.0C 

1.00 

1.09 

l.OC 

1.00 

1 

.  3  v 

1.0  0 

1.03 

1.30 

1.30 

1.0  0 

1.03 

1.03 

1.0  0 

1.00 

1.00 

1 

.00 

1.0  0 

1.00 

1.00 

1.30 

1.3  3 

l.CC 

l.CO 

l.CO 

1.33 

1.00 

1 

.03 

1.00 

1.03 

1.00 

1.03 

l.CC. 

l.C  3 

1.00 

l.CO 

1.00 

1.00 

1 

.00 

1.00 

1.03 

1.00 

1  .00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1 

.00 

1.00 

1.00 

1.0  3 

1.30 

1.0  0 

l.CC 

l.CO 

l.CO 

1.30 

1.00 

1 

.00 

1.00 

1.03 

1.00 

1.0  0 

1.00 

1.00 

1.00 

l.CO 

1.30 

1.00 

1 

.00 

1.00 

1.0  C 

1.03 

1.30 

1.03 

l.CC 

1.00 

l.PO 

l.CO 

1.00 

1 

.00 

1  .00 

1.00 

1.00 

1.30 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1 

.00 

1  .00 

1.03 

1.30 

i .:  3 

1.00 

1.00 

1.00 

1.00 

1.00 

1.30 

1 

.00 

1.03 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

l.CO 

1.00 

1 

.00 

1.00 

1.00 

1.00 

1.33 

1.00 

1.00 

1.00 

l.CO 

l.CO 

1.03 

1 

.00 

1.00 

1.00 

1.00 

1.30 

1.00 

1.00 

1.03 

l.CO 

1.00 

1.00 

1 

.00 

1.00 

1.03 

1.03 

1.33 

1.30 

1.00 

1.00 

l.CO 

1.00 

1.00 

1 

.00 

1.30 

1.30 

1.03 

1.0  3 

1.00 

1.00 

1.00 

l.CO 

l.CC 

1.03 

1 

.00 

1.00 

1.00 

1.00 

1  .  3  0 

1.00 

1.00 

l.CO 

l.CO 

1.00 

1.00 

1 

.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1 

.00 

1.30 

1 .00 

1  .0  3 

1  .  w  3 

1.00 

l.OP 

1.00 

l.CO 

1.03 

1.00 

1 

.00 

1 .30 

1.00 

1.00 

1.30 

1.00 

1.00 

1.00 

l.CO 

1.00 

1.00 

1 

.00 

1.00 

1.0  0 

1.00 

1.00 

1.0C 

1.30 

1.00 

1.00 

1.00 

1.00 

1 

.00 

1.0C 

1.00 

1.00 

1  .00 

1.00 

1.00 

1.00 

l.CO 

1.00 
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0.00 

0.00 

0.00 

0.00 

0.0  0 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

o.oc 

0.00 

0.00 

0.00 

0.00 

0.00 

10.00 

10.00 

0.00 

0.00 

3.00 

0.00 

o.oc 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

3.00 

0.00 

0.00 

0.00 

10.00 

10.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

O.GO 

0.00 

0.00 

0.00 

3.00 

0.00 

0.00 

10.00 

10.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

o.oa 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

10.00 

10.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

O.OG 

0.30 

0.00 

0.00 

0.00 

0.00 

0.03 

o.co 

0.00 

0.00 

3.00 

0.00 

10.0  0 

10.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

O.OG 

C.00 

3.00 

O.OG 

o.oc 

0.00 

C.00 

0.00 

0.00 

0.00 

0.00 

10.00 

10.0  0 

0.30 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

Q.00 

0.00 

0.00 

Ci. 03 

O.OG 

0.00 

0.00 

10.03 

10.00 

0.0  0 

0.00 

0.00 

0.00 

C.00 

0.00 

0.00 

0.00 

0.00 

0.00 

o.oc 

0.00 

0.00 

0.00 

0.00 

0.00 

o.co 

0.03 

0.00 

10.03 

10.00 

O.OG 

0.00 

0.00 

0.00 

0.00 

0.00 

o.co 

0.00 

0.00 

0.00 

0.00 

o.oc 

G.O  0 

0.00 

0.00 

0.00 

0.00 

0 . 0  0 

0.03 

10.00 

10.00 

0.00 

0.00 

0.0  0 

0.00 

0.00 

C.00 

0.00 

0.00 

0.00 

0.00 

0.30 

0 . 0  0 

0.30 

0.0  0 

0.00 

0.00 

0.00 

0.00 

0.00 

10.00 

10.00 

0.00 

0 .00 

0.03 

0.00 

0.00 

C.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

o.co 

0.3  0 

0.00 

0.00 

o.oc 

o.co 

0.03 

10.30 

10.00 

0.00 

0.00 

0.00 

0.00 

o.oc 

0.33 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0  .  0  u 

0.00 

C.0  0 

0  .GO 

o.oc 

O.OG 

10.00 

10.00 

0.30 

0.00 

0.00 

0.00 

C.00 

O.OG 

O.CO 

o.oc 

o.oc 

0.00 

0.00 

0.00 

0 . 0  L* 

0.00 

o.oc 

O.OC 

0.00 

0.00 

o.oc 

1  0  .  C  0 

10.00 

C.00 

0.00 

0.00 

0.00 

o.oc 

O.OC 

0.00 

0.00 

0.00 

0.00 

0.00 

C.00 

0.30 

0 . 0  j 

C.00 

o.co 

o.oc 

0.00 

0  .  u  0 

10.00 

10.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0  .00 

O.OG 

0.00 

o.co 

0.30 

3.00 

0.00 

C.00 

0.00 

D.CO 

0  •  0  0 

0.00 

0.00 

O.OG 

10.00 

10.00 

o.oc 

G.00 

C.0  0 

0  .  c  c 

0.00 

0.00 

0.00 

o.oc 

O.OG 

0.30 

0.00 

3.00 

0.00 

c.co 

0.00 

o.co 

0.00 

3.0  0 

0  .  0  u 

10. GO 

1  0  .  0  G 

0 . 3  0 

3.00 

0.00 

C.00 

o.oc 

0.00 

0.03 

D  .0? 

o.oc 

O.OG 

0.00 

0.0  3 

o.oc 

0.00 

0.00 

0.00 

0 . 0  u 

O.PG 

C.00 

1  j  .  0  j 

1  G  .  C  u 

O.OC 

C.00 

G.OG 

0.03 

0.00 

0.00 

C.00 

0.00 

0  .U  0 

C  .QC 

G.0  0 

o.oc 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

13.00 

1  G  .  0  C 

10. C  0 

10. GC 

10.00 

10.30 

10. CO 

10.00 

10.00 

10.00 

10.03 

10.00 

10.03 

i  o .  a  c 

10.00 

10.00 

1G.G0 

1G.00 

10.00 

1C. CO 

10.00 

10.00 
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APPENDIX  B:  FINITE  ELEMENT  PROGRAM  -  DOCUMENTATION  AND  SAMPLE  INPUT  AND  OUTPUT 
Program  FEHEAT 


1005 

C 


28 

C 


3 

1  000 
C 


4 

1001 

c 


this  PROGRAM  SOLVES  2  DIMENSIONAL  steady  state  heat  transfer  problems 
USING  THE  FINITE  ELEMENT  METHOD.  (ONLY  TRIANGULAR  ELEMENTS  MAY  BE  USED) 

DIMENSION  THE  MATRICES  WHICH  WILL  NOT  BE  STORED  IN  COMMON  STORAGE 
DIMENSION  NMJI300)*NECI300) 

DI0IMENS10NTTHE5REMAINING  MATRICES  INTO  BLANK  COMMON  BLOCK  STORAGE. 
COMMON/M1I/1MAT  <300) 

COMMON/M1R/TKM(300).GI(300) 

COMMON/M12I/NC1<300>*NC2(300) 

COMMON/M12R/X(175). YC 1 75 > » H ( 3 0 0 » . F Q I E 1 1 75) 
C0MM0N/M2I/MC1*HC2«NN*NCASE.NBHF(175) 

C0MH0N/M2R/BHF< 175) «TAI300) .RHSTC175) .RHSC300) 

COMMON/M23I/MC.NBTI 175) 

COMMON /M23R /R 1 C 175) 

C0MM0N/M3I/NE,NCN*NDF.NSZF.NBAND 
COMMON/M3R/TMI30, 175) 

C0MM0N/H13^ /NODE! 30C.3) 

COMMON /N13R/TMEI3.3) 

OPEN  THE  NECESSARY  DATA  FILES. 

CALL  C0NTRL(1,»NPDXXX*,5> 

CALL  C0NTKL<1«*E0XXXX**6) 

CALL  C3NTRL<1*MCTXXX*.7) 

CALL  CONTRL (I.'ImGXAX* .8) 

CALL  C0NTRL<1,*SHFaXx**9) 

CALL  C0NTRL<l.*ChSXXX**14) 

CALL  CONTRL  <1 « *uUANXX*tll) 

CALL  C3NTRL<2. ’OUTPUT*. 12) 

CALL  C0NTRL(1»’TI0TXX**13) 

READ<11»107.ENO=10051NN.NE.MC«MJC*MC1.MC2.NIT 
COMT I NUE 

-  Ik -  all  MATRICES  ANO  PARAMETERS  TO 


ZERO, 


INITILIZE 
NCN  =  3 
NOF=  I 
NCA$E=1 
NSZ  F  =NN  *NDF 
TKM<1)=.833 
DO  2  I=l,NN.l 
X(I )=0.0 
YCI)=0.0 
NBT(I)=0 
BT (  1  )  =  0 . 0 
NBHFU  )  =  0 
BHF ( I ) =0 . 0 
RHS ( I ) =0 . 0 
FQIE(I)=0.0 
CONTINUE 
DO  28  U=l.NE.l 
NCI <  U  >  =  0 
NC2I J)=0 
HC J)=0.0 
TA(U)=0.0 
IMAT<J>=0 
NMJ< J)=0 
Q I <J)=0.0 
CONTINUE 

READ  IN  THE  NODAL  POINT  DATA. 

REA0(5*101*ENO=1000>  I»X<II»V(l> 

CONTINUE 

CONTINUE 

READ  IN  THE  ELEMENT  DATA. 

READ<6*iS2*END=1001)  J»N0DEIJ*l)»N0DEtU»2).N0DEtJ*3>*IMATIJ> 
CONTINUE 

Tcl'oT*  THE  FIXED  TEMPERATURE  BOUNDARY  CONDITIONS. 
IFIMC.EO.O)  GO  TO  5 


U^OOwi  ►—(i-g  O'f*)' 


CO  h  1*1=1. KC.l 

REALM  7, 1  C3»  END=1C  02  )  M  .  N  ti  T  (  N )  » tiT  (  N  HT  (H)  > 
k  LOMHUl 

1002  CONT I NUE 

r-.i-AO  IM  THE  INTERNAL  HEAT  GENERATION  VALUES  FOR  EACH  ELEMENT. 
1FIMJC.EO.0)  CO  TO  1 
'JO  7  I*J  =  1,IRJC«1 

RLt  Mb,  105)  MJ,NMJ(MJ),QI<NMJCMU>> 

UMIVJt 

KE AD  I M  THE  HEAT  FLUX  VALUES  FOR  NODES  WHERE  HEAT  FLUX  IS  SPECIFIED 
IF  (MCI.  Ecu)  CO  TO  b 
1)3  ■*  IMC1=1,MC1,1 

KEaMS.ICu)  «l.N:iHF  (Ml  >  ,UHF  (NbHF  <M1>  ) 

CONTINUE 

cOh  ‘-OuNJARY  SEGMENT.;,  SUBJECT  TO  CONVECTIVE  HEAT  TRANSFER  READ 
IN  COaRESPJNJJ N3  NODES  »  h  » A  N  0  AMblENT  TEMPERATURE. 

IF  (MCi'.t  ...0  )  CO  TO  1C 
OO  11  M2=i.MC2,l 

K  r  (  l  4 , 1  3  1.  >  NE  C  C  -I?  >  .  ,\IC  1  (  NEC  t  v?  I  )  ,  NC2C  VEC<  M2  )  )  ♦  H<  NEC  I M?  >  >  . 

a  i  a  r-i  e  c  (  m  ?  n 

11  utJNTINJF 

1  '  C ON T I SUE 

*  •■  1  I  r  (  1  i  »  1 «'  C  >  N".*Nf' 

.  i  t  f  1 1 ; , l  ?  n 

0  3  '-l  1=1,\N,1 
.  "  I  r  t  <  1  „  ,  1  M  )  I  ,  \  (  I  )  ,  Y  (  I  ) 
il  CONTINUE 

.RITE (12 .122) 

JO  H  2  J=l»Nt.l 

N M  T  L ( 1 3 , 1 L2 ) J ,  .0 0:M U. 1 > ,NG3C  <  J, ? ) .NODE  <  J. 5) » I M  A T ( J> 

( 2  '  OMlNJf 

I  F  (  MC  )  ■)  i  «  y  1  .  3  1 

7  1  *HlTf(l2.12S) 

uO  a 5  ^ -  1 « f’ 2. 1 

»'F  1  TL  (  U  .  1  C  i)  !-,  vO  (  ")  ,  T(N'  T  <  M>  > 
l1:  3  0  N  T  1  'I  J 

72  LOnTINJE 

IF ( HJC ) 34  •  >4, 32 

7;  .  k  i  t  r  ( i  .  i  r  j  ) 

Jo  --4  I=1,MJC,1 

4*1  TE(  12  .1  33>MU.N«J(>IJ»  *01  (NMJ(MJ)  ) 

•)A  CONTINJE 

J1*  CONTINUE 

IM  MCI  >'■:*.  ,s,  -3  7 
•7  7  .KITE  112 .120) 

3  0  "j  ■11  =  l,*'CI.l 

w  k  I  T  C  <  1  2  *  1  0  3  >  M 1  .  N :  •  H  f  (  *  1  )  .  T  H  F  »  N ; » H  F  (  M  1  )  ) 
u'  CONTI N  JF 

7)  CONTINUE 

IF  <'  C  2  )  ~>c‘ 

7  -s  .  i\  1 1  f  ( u , i ;  t- ) 

UO  -in  *2=1, PC?, I 

Wn]  TE(l^»12h)M2(w2)»NCHNtC(M?>>.NC2CNtC(M2))«H(NECCM2)), 

a  T  M  NIC  <  M?>> 

UK  COM  I  NJf 

7  k  r  o  n  t i n  j  e 

C  FINU  K4N0.I3T-  FT,.  TkIANGJLAR  EILWENTS. 

!’G  tb  Isl.HE.l 

Ll  =  I  AlS<N0..E(  I.1>-N0DE<1.3>) 

L  2  =  1  At-  0  (NOuf  <  I  ,  j)  -\  CUE  (  I  .2  )  > 

L  J  =  I  At  S<  NO  OF  <  1  .  T-NGDE  I  I  .1  )  > 

Na=PAX  i  (  MX  »L  1  »L2  »L  .'  > 

,  b  CONTINUE 

M  r  i  N  J  =  'U'  F  »  <  M  X  ♦  1  ) 

4  a  IT F (12.12") N 'AN 3 
.■•U."  =  M  A  s!>-  I 

i  J-: 

1  0r  =  i 

CALL  FO*** 

CALI  Ff.HM  I) 

CALL  MC-tJ<i\bS»T',«NS2F»l*MUU»IOP.O»E'E“<>»lER»IQ) 

W* I TE( 12.127) 

JG  200  1=1. V-i, l 
.kIIl(IL.IJL)  I  ♦  K  Ft  S  (  I  ) 

200  CONTINUE 

IF(MT.CC.J)  GO  TO  22  8 

CALL  IS0THV(8H'>. NODE. X.Y.NE. NIT) 

22k  CONTINUE 

1 j 1  FORMAT ( 1 1>  «  2F  1  0  .  A  ) 

132  FORMAT  ( lb, <*16) 

1  J3  FORMAT* IN, 16, F1C. A  ) 

l  06  FORMAT  (JI6.2F10.4) 

137  FOnMAT ( 7Ib> 


123 

129 

125 

126 
127 

130 


C 

C 

c 


3*0 


X 

X 

X 

X 

X 

X 

X 

X 

X 


F0RMAT(5X»'THE  BANDWIDTH  IS', 16) 

FOR  M AT ( /  / »5X  » ’TOTAL  NUMBER  OF  NODES  ='.I6,3X» 

♦TOTAL  NUMBER  OF  ELEMENTS  =  •»16> 

FORMAT (//  ,2X, *NODE',3X, 'GLOBAL », AX. 'GLOBAL*  ,/»3X,'N0.» ,6X,'X», 

9X,  »Y  *  ,1  X  »/,  • - -  » 

FORMAT ( // ,1 X, 'ELEMENT • *  1 X , 'NODE » »2X , 'NODE '» 2X , 'NODE • ,2 X » 

JMATL. »,/,3X, 'NO.', AX, '1 ».5X,'2',5X,'3» , A X, »T YPE • ,/ , 1 X* • - ♦» 

FORMAT C//.3X. 'B  .C  .  •  »2X  «  •  NODE  •  *  2 X  ,  •  TE HP •  , / »  A X, 'NO.'.SXt  »NO.*»2X» 

•(F)',/, IX,' - •) 

FORMAT (//  2X, »S.C. • »2X • 'NODE • * 6X « »Q I • ./ ,2X • 'NO. • «3X * 'NO. * «5X* 
•(UNITS)') 

FORMAT  (  // ,2 X , 'U«C'«2X,'N0DE**6X**HF**/,2X,*N0.'»3X**N0«*«5X, 
'(UNITS)') 

FORMAT  (  // »  1  X, 'SEGMENT '.IX,  'NODE* ,2 X , 'NODE • » AX , • H • , 3X , • AM6 » ,/ , 3X, 
•NO •'•AX* *1 •  «  5  X  , *2'«3Xf' (UNITS) •  .2X, 'TEMP(F)  •  ) 

FOR M AT (//« 1 3X .' NODE • ,2 X.' TEMPERATURE  (F)», 

/. 12X «  * - •) 

FORMAT ( 1 0  X  , I6.F13.2) 

CALL  CONTRL ( A , • NPDX XX • « 5 ) 

CALL  CONTRL (A,'EDXXXX'»6) 

CALL  CONTRL (A »'UCTXaX'*7) 

CALL  CONTRL ( A • • 1HGXXX • .8) 

CALL  CONTRL ( A , »SHF X XX • *  9 > 

CALL  CONTRL (A ,»CBSXXX • » 14} 

CALL  CONTRL (A.  "UUANXX',11) 

CALL  C0NTKL(4, 'OUTPUT', 12) 

CALL  CONTRL (A, 'TIOTXX',13) 

CALL  EXIT 
END 


SUBROUTINE  rSM(R) 

0  I  MENS  I  ON  XC( 3) , YC( 3), T< (300) 

COMMON/M1I/IMAT(30C> 

COMMON/M1R/TKM(300)»QI  (300) 

COMMON /Ml 2  I /NCI ( 3 0 0 > , NC 2 ( 3 0 0 > 
COKMON/M12R/X(175)»Y(175)«H(30O)»F3IE(175) 

C0MM0N/M13I/N00C( 300,3) 

COMM3N/M13R/TME(3,3) 

<  IS  THE  ELEMENT  NJMuER. 

T<(<)=TKM( IMAT(K)  ) 

Nl  =  .NOOE(K.l  > 

N2=  NOOl (K,2  > 

N3=N00E(<.3 > 

DEFINE  THE  ELEMENT  NOOAL  X  A‘,0  Y  GLOBAL  COORDINATES. 

XC(1)=X(N1) 

YC(1)=Y(N1) 

XC<2)=X(N2» 

YC(2)=Y(N?) 

XC(3)-X(N3) 

YC( J)=Y(N3> 

A  =  1  .0 
A  A  =  1 , 0 

DEFINE  THE  A'S,ft'S,ANO  C'S  OF  THE  LINEAR  INTERPOLATION  FUNCTIONS. 
Dl=YC(2)-YC(3) 

B2=YC(3)-YC(1 > 

B3=YC(1)-YC(2) 

C1=XC(3)-XC(?> 

C2=XC(1)-XC(3> 

C3-XC ( 2 ) -XC  ( 1 ) 

DETERMINE  THL  ELEMENT  AREA 

3EL  =  ABu(0.5»(XC(l)*( YC(2)-YC( 3 ) ) 'X C ( 2) • ( YC 1 3 )-YC 1 1 > ) *X C I  3) 

X  *(YC(1)-YC(2) ) ) ) 

FORM  THE  INFLUENCE  MATRIX  FOR  TEMPERATURE. 
CONST=(TK(K)*A/(A.O*OEL))*AA 
TMEt  1,1  >  =  (bl*»2.Cl"2  >  'CONST 
TNE<1»2)=(B1*0?'C1*C2)*C0NST 
TME(l»3)  =  (‘Jl'l‘3'Cl*C3)*CONST 
TME(2»1)=TME(1,2) 

TNE(2,2>=<B2**2'C2"2>  'CONST 
TME(2»3)=<32*B3'C2«C3)*C0NST 
TME(3»1)=TME(1,J) 

THE (3.2)=TMt (2,3) 

THE ( 3,3 >=<B3»»2'CJ»*2> .CONST 
IF(NCl(K))3lb,330,3AO 
KK1=NC1(<  ) 

KK2=NC2(<) 

IF(KKl.EO.Nl)  <1  =  1 
1 F ( KK1 . E  0 . N 2 )  <1  =  2 
lFiKKl.EG.N3)  <1  =  3 
1 F ( <K2  »£G  »N1 )  <2  =  1 
IF (KK2.EQ.N2)  <2=2 
IF(KK2.EG.N3)  K  2  =  3 


71 


330 

c 


315 

320 

10 


C 


780 


V- 

7S0 


SCO 

•110 

c 

8  2  0 


S3? 
8  4C 


700 

c 

c 

c 


c 

c 

c 


300 

c 

c 


c 

c 

c 


IF(<i.£Q.K2»  SO  TO  315 

CONSTC  =  H<K>*<$ORT  C ( X < N Cl (K ) > -X (NC2 (K > > > * *2* t Y ( NCI (K)  ) - 
X  V  <  NC2  (K  > ) ) **2) > 

TME(Kl»Kl)=TMECKl,Kl)»C0NSTC/3. 

TME(Kl»K2>=TM£(Kl,K2>*CONSTC/6. 

THE(K2«K2>=TH£tK2»K2)*CONSTC/3. 

TMEtK2,Kl>=TME(K2.Kl>*CONSTC/6. 

continue; 

FORM  INFLUENCE  MATRIX  FOR  INTERNAL  HEAT  GENERATION. 
QIt=QI <K> 

C0NSTQ=0EL/12. 

QN=4.*QIE»C0NSTU 
FQIE<N1>=F3IE(N1>*6N 
FJIE(N2)=F31E(N2)*QN 
FOIE  (N3)=FaiE(N3)*QN 
GO  TO  12 

WR  t  TE ( 12, 320 ) 

FORMAT< 20 X#»»***»« ‘ERROR ***»**•) 

CONTINUE 

RETURN 

END 


SUBROUTINE  FRHStBT) 

DIMENSION  t)L  (  3  C  0  )  »  R  :'L  (  3  0  0  ) 

DIMENSION  8 T < 1 75 ) 

C0KM0N/M12I /NCI  1303 > ,NC2 I  30 0 > 

C0MM0N/M12R/XU75).  Y C  1  T5 >  . H C  3  C  0  >  «F  U 1 E <  1  T5> 
C0MM0N/M2I/HC 1,MC?»NN»NCASF,NF HF( 1 75  > 

C0M»0N/M2R/BHF(  175)  ,TA(  300  > . R HS T ( I  75  I , SHS « 3 j 0 ) 
C0MM0N/M231  /MC,NbT (175) 

CUMM0N/M23K/RH  17h> 

DO  7  SO  1=1«NN,1 
RHS ( I ) -F  3 1 C ( I ) 

CONTINUE 

IF(MC1)h10»810,790 

INSERT  THE  BOUNDARY  hEAT  FLUX. 

DO  800  1  =  1, MCI, 1 

HMS(NcHF (I > )  =  KHS(NbHF  U ) J-HHF  (NhHF  (  I  )  ) 

CONTINUE 

IF (KC2)8AL,h40,’2C 

ACCOUNT  ^OR  CONVECTION  AT  THE  BOUNDARY . 

DO  S3C  1=1, R C 2,1 

bL< I)=SuRT< (X (NCI ( I ) )-X(NC?( I )) >**2MY(NC1( I > >-Y(NC2(I 
K8L(]>=2.»3.1Alj9*<(0.5*(X(NCl(I)>*MNC?U)>>>l**? 

I  F  (  NC  A  SE  »E  >«  •  2  >  BL(I)=Ki>L(l)*“L(I> 

Cl  I NF  =  H( I )»T A(I >*nL( I ) /2.0 
RHS(NC1 ( I >  >=RHS(NC1 (I) )*CTTNF 
Kh.sT  (NL1(  I  )  )  =  RHoI  (NC1(  I )  )*CT  1  ‘.p 
RHS(NC2(D)=RHS(NCr<I))»CTINF 
RHST  (NC2<  m  =  RHST(NC2(I  )  >  ♦C  T  I  NF 
CONTINUE 
CONTINUE 

INSERT  T  H  t  BOUNDARY  CONDITIONS  ON  TEMPERATURE  • 

DO  8  0  0  ';  =  1,"C,1 
I  =  N  8  T  <  \  ) 

KHS  (  1  )  =K  1  (  1  >*I.T  (  I) 

CONTINUE 

RETURN 

END 


SUBROUTINE  FORM* 

FOR'-’S  STIFFNESS  MATRIX  IN 

UPP*  R  TRIANGULAR  F  O  M 
COMM ON/M231/MC.NB 1(1751 
C  0  M  M  0  N  /  M  2  3  R  /  K  1(175) 

COMMON /M .3 1  /NE  ,NCN  ,N!)F  »NS  2F  ,N8  AMi 
C0MM0N/M3K/TM ( 3  J  » 1 7  j) 

COMMON/M  13  I /NOin  ( 300,3) 
C0VM0N/M13K/TMF;(3»3) 

DIMENSION  ST ( 52  sO  ) 

EOU1 VALENCE  ( ST (  1 )  ,T M (  1  , U  ) 


ZERO  STIFFNESS  MATRIX 


DO  JOG  N=i,NS2F 
UO  3  CO  «=l,NfSANU 
TM(.M,.N)=0 

SCAN  ELCMENTS 


DO  A10  N  =  1 , N C 
CALL  TSM(N> 


RETURN  ESTIFM  AS  STIFFNESS  MATRIX 
STORE  ESTIFM  IN  SK 
FIRST  RURS 


)  >**2> 


00  350  JU=1.NCN 
NR0W8  =(NODE<N, UJ)-1)*NDF 
DO  350  J=1,N0F 
IF  (NROWB)  350,305,305 
305  NROWB  =NR0WB  *1 

I=« JJ-1>*NDF*J 

C  THEN  COLUMNS 

C 

DO  330  KK=1,NCN 
NC0LB= ( NODE <N ,KK > -1 > *NDF 
DO  320  K  =1 »NDF 
L  =  I KK- 1 >  *NDF *K 
NC0L=NC0LE*K*1- NROWB 

C  SKIP  STORING  IF  BELOW  BAND 

C 

C 

1FINC0L)  320, 323, 3U 

310  TMINCOL. NROWB)  =  TH<  NCOL, NR  OWB )  ♦  TKEFI.L) 

320  CONTINUE 

330  CONTINUE 

350  CONTINUE 

400  CONTINUE 

C  INSERT  BOJNDART  CONDITIONS 

C 

DO  500  N=1,MC 
I=N8T(N> 

TMU,I  J=TMU,I)*1.E15 
RU1):TH(1,  I) 

500  CONTINUE 

DO  1  J=1.NSZF 
169  FORMAT (’ROW  *,I2> 

166  FORMAMIOE  10.21 
167  FORMATE//) 

1  CONTINUE 

C 

C  CHANGE  STIFFNESS  MATRIX  TO  CONSECUTIVE 

C  STORAGE  LOCATIONS 


M  =1 

DO  777  I=1,NSZF 
K  =NSZF-I  *1 
00  700  J  -1 ,NDANQ 
IF  CJ-K)  705.705,777 
735  STEM)  =TM(J,1) 

700  M  rrt  ,1 

777  CONTINUE 

198  FORMAT  I  16. 6X  ,  It) 

RETURN 

END 

C 

C 

c  - 

SUoROUT INC  KCHB ( R ,A,M,N»MU3, I OP.EPS.IER, 101 
C  *  I  SUBROUTINE  MCR9  <  SE  E  IBM  SSP  ROOK).  USED  TO  SOLVE  MATRIX 

C  EQUATION  WHEN  5QU ARE  MATRIX  IS  S TMME TR I C , BANDED, AND  POSITIVE 

C  DEFINITE.  NO  NEW  INPUT  DATA  REQUIRED. 

DIMENSION  R C300),A(525C) 

DOUBLE  PRECISION  TOL.SUM.PIV 
IF<IA3S<I0P>-3>  1.1,43 

1  I  F  E H JD )  45,2,2 

2  NC=MUO*l 
IFCM-MC»46,3,3 

5  MRsM-MUD 

I  ER  =  0 

I F  E I  OP >24, 4, 4 

4  IEND=0 
LL0ST=MUD 
DO  23  K  =  1.M 
I  ST  =  I  END* 1 
IENO=IST*MUD 
J=K-MR 

IF (J)6,6,5 

5  I END= 1 END-U 

6  IFCJ-118,8,7 

7  LLDST  =LLDST -1 

8  LMAX=MUO 
U=MC-K 

IF(J)10,10,9 

9  LMAX=LMAX-J 

10  1D=0 
T0L=AE1ST>»EPS 
DO  23  I  =  I S T  ,  I  END 
SUNrO.OO 
IF(LMAX>14»14»11 


11 


12 

13 

14 

15 

16 
17 
16 
19 


20 

21 

22 

23 

24 

25 


2b 


27 


2  6 
2  9 

30 

31 

32 

53 

34 

55 

36 

37 


38 


39 

40 

41 

42 

43 


LL  = 1ST 
LLD=LLOST 
00  13  1.  =  1  »L  MA  X 
LL=LL-LLD 
LLL=LL*10 

SUM=SUH+A(LL) *A(LLL) 

B  =  SUH 
B=ABS ( B ) 

IF(B.LT.1.0E-35>  SUM  =  0.000 
IF(LLD-MUD>12,13.13 
LLD=LLD*1 
CONTINUE 

SUM=OBLE(A(I ) ) -SUM 
IF(I-IST>15,15.20 
IF(SUM)47»47,16 
IF (SUM-TOL  >17,17,19 
IF(IER)18,18.19 
IER=K-1 

P I X  =  DSQR  T  < SUM  ) 

A  < I ) =P 1 V 
PIV=1.0U0/PIV 
GO  TO  21 

AtI)=SUM«PIV 
I  0  =  1 0*1 

IF< I0-U>23,23.22 
LMAX=LMAX-1 
CONTINUE 
IF ( IOP  >  24,44,24 
I U=N*M 

IEN0=IABS< IOPI-2 
IF( IEN0I25, 35,25 
I S  T  =  1 
LMAX=0 
U  =  -NR 
LLJ$T="UO 
JO  54  K  =1  ,  M 
P I  V  =  A  < 1ST) 

IF<PIV>26. 4fa»26 

PIV=1.0O0/PIV 
00  30  I=K,ID,M 
SU8=0.000 
IFtLMAX)30,3C,27 
LL=I ST 
L  LL  =  I 
LLO=LLDST 
uO  29  L=1,LMAX 
LL=LL-LLO 
LLL=LLL-1 

SUM  =  SU>1*A<LL)  *R  <LLL  ) 
IF(LLD-MUD)2R, 2-»,29 
LL0=LLU»1 
CONTINUE 

R<I>=PIV*10»LE(M I) >-SUM) 
IFC'5C-X)32.32,31 
LMAX=K 
I  ST  =  I S  T  ♦  HC 
J  =  J»1 

IMJ>34, 34,33 
1  ST  =  I S  T - J 
LLJST=LID3T-1 
CONTINUE 

IF 1 1 EN 0 ) 35 , 35 , 44 

lST  =  M*<MU0«<M*M-f'Cn/2*l 
LMAX=Q 
K  =  M 


I ENO= 1ST -1 
ISTrlENO-LMAX 
P I  V  =  A  <  1ST) 
IF(PIV)37,49,37 
PIV=1.0OQ/PI 


V 


L  =  I ST+ 1 

00  40  I =K , I  0, M 
SUM  =  O.OOCi 
IF<LMAX>40, 40,38 
LLL  =  I 

DO  39  LL=L , IENO 
LLL=LLL*1 . 

SUM=SUM*ACUL)*R (LLL I 
R  ( I  >  =P  I V  » (  OBL.E  ( R  (  ID-SUM) 
IFIK-MK)  42,42,41 
LNAX=LMAX*1 
K  sK  *  1 

IF<K>44, 44,36 
IER=-1 


T 


m 
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GO  TO  44 

♦  5  ltR=-2 

GO  TO  44 

46  IER=-3 

GO  TO  44 

47  IER=-<K*10> 

WHITE  110,9000)  SOM*  A  ( 1 ) 

9000  F0KMAT<5X»C14.7»E14.7> 

GO  TO  44 

48  I ER  =  -5 

GO  TO  44 

49  I£R=-6 

44  RETJRN 

END 

C  - 

SUBROUTINE  ISOlHM<T.NODE»XC*YC,NE,NIT) 

01  HENS  I  ON  MODCt ?0  0,3>»T(300),XCU75>,YCU75> 

DIMENSION  KK<13>,TISO<1O),X<10O>,Y<1QO>»XS<1OO>,VSI1OO> 

00  350  I  =  t«.\IT«l 

J=0 

HEAJI13«230  .ENJ  =  190  ))  TISOU) 

230  FORMATS. 2) 

TT=T  1S0(  1 ) 

DO  57  *=1,NE,1 
Nl  =  NOOE  <<  « 1  ) 

N2  =  N00E  » 

N5=NODE<R,31 

Tl-TTNl) 

T2-TIN2) 

T  3-T (\ 3  ) 

Xl=XC(Nl  ) 

Y  U  T  C I  N  1  ) 

X  2  =  X  C  <  4  J  ) 

Y2=YC<N? J 
X3  =  *C<  N  3 ) 

Y3=YC< x 3) 

701=11-12 
TJ2=T1-T J 
T03=T2-T3 
lF(TDl)H,6i),o 

60  IF(A?S(T1-TT)  .oT..311  GO  TO  ?  • 

J  =  J«l 

X<J)  =X1 
Y< J)=Yl 
J  =  J»1 
X  (  J  )  =  X  ? 

y  «  j  >  =  y  ? 

«  0  TO  £ 

t  t  h  i  =  t  ;■ 

T  C  L  =  T  1 

IF<TT-TCl>2fc,13,10 

10  IFITH1-TT)2V«11,11 

11  J  =  U*1 
K4TT0=<TH1-TT»/UH1-TC1» 

XIJ)=X?-(<x2-X1)»4ATI0> 

Y<J)=Y2-MY;'-m«RATI3l 
GO  TO 

9  I41=T1 

TCI =T2 

lF(TT-TCi>3>-.,lv,i2 

12  1FITHI-TT)2>,13,13 

13  J  =  J*1 
RATI0=<TH1-TT  I/U41-TC1 1 
x(J)=Xl-<<Al-X2>»RAT10> 

Y<J)=Y1-HY1-Y2>*RATI0> 

28  IFU02)  14,61 ,16 

61  IFt*tJSm-lT>.GT..ei>  oc  TO  20 
U  =  J*1 

X  <U)=X 1 

Y(J>=Y1 
J  =  J«1 
X< J)=X3 
YI J»=Y3 
GO  TO  23 

14  TH2=T3 
TC2=T1 

IF( TT-TC2129,  16,16 

16  IFITM2«TT>29,17,17 

17  J  =  J»1 
RATI0=(TM2-TT )/<TH?-TC2) 

X<J)=X3-<(X3-X1>*RATI0> 

Y(J»=Y3-UY3-Y1  »««UTIO) 

GO  TO  29 

15  TH2=T1 
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TC2=T3 

IF(TT-TC2>29»18»18 

IF(TH2-TT>29,19,19 

U=J*1 

kATI0=(TH2-TT  >/<TH?-TC2> 
X<J)=X1-UX1-X3>*RAII0> 

Y(J>=Y1-<(Y1-Y3)*RATI0> 

IF(TD3)20«62«21 

IFCAE>S(T2-TT> .GT..01)  GO  TO  30 
J  =  J*1 
X<  JJ=X2 
Y(  J»=Y2 
J  =  J*1 
X<  J)=X3 
Y(J)=Y3 
GO  TO  30 
TH3=T3 
TC3  =  T2 

IF(TT“TC3>30»22»22 

1FITH3-TT) 30,23,23 
J=J*1 

RATI0=<TH3-TT>/CTH3-TC3> 

X(J>=X3-((X3-X2)«RATI01 

Y<J)=Y3-<(Y3-Y2)*RATI0) 

GO  TO  30 
TH3=T2 
TC3=T3 

IF(TT-TC3)30,2A,2A 

IF<  TH3-TT>30«25,C5 
J  =  J*1 

RATI0=(TH3-IT)/(TH3-TC3) 

X(J>=X2-< <X2-X3I*RAT 10) 
Y(J)=Y2-((Y2-Y3)*RATI0) 

CONTINUE 

CONTINUE 

K  =  1 

XSl 1 )=X ( 1 ) 

Yi( 1)=Y( 1) 

00  PO  M1=2,J,1 
ki-: 

DO  *'?  M2=l,K,l 

IF  t  At*S(  Xs,< '*2> -X  (Ml)  >  .GT  .  .0  II  GO  TO  V  0 
IF<«uS(YS(M2>-Y(Ml>>.GT..01>  GO  TO  TO 
K  I  -  X  I  ♦  1 

CONTINUE 

IflKI.If.SI  GO  TO  r.  i 
K-K*l 

XS(X>=X(M1> 

YS(X)=Y(M1  ) 

FONT INUT 
K  K  (  J  )  —  < 

mRITE(12«210)TT 

FORMAT!//, X, ♦COORDINATES  FOR  THE*,F8.2,t 

X  - - ♦) 

KK  <  1  )=K 

OG  32b  N-l »R  ,  1 

WRITE  (  12,220  XS<N),YS(N) 

CONTINUE 

FORMAT (1CX,2F1Q.A) 

CONTINUE 
CO.NT  I  NUE 
RETURN 
END 


ISOTHERM*,/, AX, 
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Table  B1.  Sample  output  from  FEHEAT  —  provides  initial  conditions  and  final  temperature  distribution. 
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1.0000 

G.G003 

69 

J  .9239 

0.3827 

TOTAL  NUMBER 

;  OF  ELEMENTS  =  166 

73 

71 

0.7071 

C  .  36  27 

0.7071 

0.9239 

TOTAL  NJMoE9  OF 

NODES  =  104 

72 

73 

3.0030 

0.4142 

1.0000 

-1.0000 

74 

0.7071 

-1.000G 

75 

1.0000 

-1.5000 

GLOBAL 

7b 

1.0000 

-0.7071 

NODE 

GLOBAL 

77 

1.0000 

-Q.4142 

NO. 

X 

Y 

7ft 

1 .0  300 
l.oooe 

0.414? 

0.7071 

1 

: .  o  c  o  o 

-0. 1  OOC 

90 

1.0000 

1.0000 

2 

L  .  0  3  b  3 

-0.0924 

91 

0.7071 

1.0000 

5 

0.0707 

-0.0707 

82 

0.4142 

1.0000 

4 

0.0924 

-u . 1383 

3  3 

J.GOC0 

-1.4C0C 

c 

0  .  1  'J  0  0 

0.0000 

34 

3.4300 

-1.4000 

b 

0.0924 

u  •  C<  3  b  *5 

85 

3.0000 

-1.4030 

7 

i>  •  j  7  0  7 

0.0737 

3  b 

1.4300 

-1.4000 

ft 

0  •  0  3  b  3 

0.0924 

97 

1.4300 

-1.3C0G 

rf 

C.C300 

0.1000 

96 

1.4000 

-0.6G0C 

1  0 

0  •  o  0  C  0 

-0.1500 

89 

1*4030 

-  0 . 2  C  0  3 

1 1 

0 . 0  -i  74 

-0.1386 

9  0 

1.4300 

0.2005 

12 

-0.1061 

91 

1 .4000 

0  •  6 1  G  0 

13 

0 .1386 

-3.C574 

9  2 

1.4000 

1  •  G  0  0  0 

14 

l  •  1  D  i<  J 

C  •  i;  G  0  L 

93 

1.0300 

-  2  •  0  G  3  j 

1 1> 

o.i :sb 

0.0574 

94 

3.4300 

-2.0000 

lb 

0 . 1  0  fa  1 

G  •  1 G  b  1 

95 

J  .  o  3  0  0 

-2.G0CC 

17 

C  •  0  b  7  4 

0.1386 

9b 

1.4300 

—  2  •  j  G  0  U 

1  8 

•j  •  G  0  J  0 

0.1500 

97 

2.0000 

-2.0003 

1  9 

0.0100 

-0.2000 

9  8 

3 . 0  0  C  0 

-1.4000 

2  0 

0.0765 

-0.1848 

9  4 

2.0000 

-1.0003 

21 

0.1414 

-0.1414 

103 

J  e  u  j  C  G 

-0.6000 

2? 

0  •  1  c  4  6 

-G»l 7  b b 

131 

2  •  0  j  <J  j 

—  J  •  J  L1  l‘ 

23 

0.2000 

0.UG33 

102 

*  •  0  j  3  3 

D  .  2  0  0  0 

24 

0.1946 

0.0769 

103 

2  •  u  o  l  G 

C  •  t  0  0  u 

2  b 

0.1414 

0.1414 

104 

l  •  0  C  G  j 

1 . 0  0  0  0 

2b 

J  •  0  7  6 

C  .  1  8  4  4 

27 

C.0300 

0.2300 

2b 

j  •  0  j  j  0 

-  0 . 3  C  C  J 

ELEMENT 

\  O  J  £  \UL-£  VODL 

M  A  T  L  » 

2  < 

t  .1 14fc 

-0.2772 

NO. 

i 

3 

TYPE 

-  0  «  2  1  2 1 

D  J 

3  i 

0.2772 

-0.114b 

1 

10  3 

1 

1 

3  7 

0.3000 

0  e  «i  u  J  l 

10  11 

2 

1 

33 

1.2772 

0 .  11 4  8 

3 

1 1  1 

-> 

1 

34 

0.2121 

o  .  i  2  i 

11  1.: 

3 

1 

3d 

0.1144 

0.  ‘7  7? 

12  4 

3 

1 

3  fa 

0.0003 

c.  3  0  1  C 

6 

12  13 

4 

1 

37 

J  *  (i  J  wi  J 

-0.4200 

7 

13  5 

4 

1 

3  b 

..1607 

-  0  e  3  r.  r 

8 

13  14 

5 

1 

39 

0.2970 

-  C  .  2  9  7  0 

9 

14  b 

D 

1 

4  a 

)  •  5  o  ft  3 

-  u .  i  fa  a  7 

13 

14  1  J 

6 

1 

4  1 

..  *4200 

3.0CC0 

11 

15  7 

6 

1 

42 

1 . 3  4  fa  0 

0.1607 

1  2 

15  lb 

7 

i 

43 

-.2470 

3.2970 

13 

lb  d 

7 

i 

44 

-.1607 

l  •  1>£  H  41 

1  4 

16  17 

6 

i 

4  b 

•.  •  U  0  «•  0 

C  .  4  2  0  0 

1  5 

17  9 

8 

i 

4b 

j  •  G  0  0  G 

-  0 • 58  0  C 

1  6 

17  1  r 

i 

47 

G  •  2  2.  i?  J 

-0.5359 

1  7 

19  11 

10 

i 

48 

i  .41  11 

-0.4101 

18 

19  23 

1  1 

i 

4b 

L .5369 

-0.0220 

19 

:o  12 

1  1 

i 

5v 

J  •  h  t  0  C 

0  •  l  0  0 1 

20 

23  21 

12 

i 

51 

C .5359 

0.2220 

21 

21  13 

12 

i 

52 

1.4101 

0.4101 

22 

21  2  2 

13 

i 

53 

- • 222  C 

0.9359 

23 

2  2  14 

13 

i 

54 

c.  •  c  c  g  c 

0.5800 

24 

2  2  2  3 

14 

i 

55 

C  .0000 

-0.7600 

25 

2  3  15 

14 

i 

5b 

0.2985 

-0.7206 

2b 

23  24 

i  b 

i 

57 

0.5515 

-0.8515 

2  7 

24  16 

15 

i 

5  fa 

0.7206 

-0.2985 

28 

24  25 

16 

i 

55 

0.7800 

0.0000 

29 

25  17 

16 

i 

6  0 

0.7206 

0.2965 

30 

25  2b 

17 

2 

61 

0.5515 

0.5515 

31 

2b  18 

17 

1 

62 

0.248b 

0.7206 

32 

26  27 

18 

1 

63 

J  •  0  0  C  0 

0.7603 

33 

2b  2  0 

19 

1 

64 

0.0000 

-1.0000 

34 

2a  29 

?0 

1 

65 

0.3927 

-  0 . 1=1 2  39 

35 

29  21 

2  0 

1 

66 

0.7071 

-0.7071 

3b 

29  !  0 

21 

1 

67 

0.9239 

-0.3827 

37 

3C  22 

21 

1 

38 

30  31 

22 

1 

39 

31 

4  :■ 

31 

4  1 

3  2 

4 

•  3 

44 

<  3 

4b 

3  4 

4  D 

34 

4  7 

3  3 

4  6 

1  2, 

4  1 

7  7 

5  0 

•  1 

5  2 

7  8 

5" 

I(  -4 

5  3 

<  r, 

5  '* 

V  9 

jb 

4 

">b 

4  r. 

5  1 

<♦  i 

5fi 

3  4 

4  r 
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4  S 
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4  4 
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4  4 

o 
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4  b 

2  7 

47 
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4  / 
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T* 
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72 

4 
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•*  it 

7  J 

..  1 

74 

.  . 
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>1 

7  6 

1 

77 

*' 

7  b 
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79 

M  r' 
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9  7 
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3/ 

!•  t> 
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.  t 
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7 

3 
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w> 
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O 
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1 
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3  ' 

t 
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4  7 
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i  ••• 
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1  D 
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7 
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V' 
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n 

22  2 

71 

1  2  2 

7  i 

1  2  J 

f  4 

224 

IS 
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D 

116 

7  7 

1  1  7 

bb 

llr 

7  r 

119 

7a 
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4  2 
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73 
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74 

23 

32 

32 

7  3 

?4 
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33 
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4  4 

34 

2  b 

2e 

;  b 

3  6 

<.  fa 

27 

,  fa 

3  L 

7  7 

2  f 

7  fa 

74 

2  9 

3  J 

29 

3  9 

30 

31 

7  3 

4  & 

31 

3  7 

3 1 
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.1 2 

3  3 
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34 

3 
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5  > 
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3- 

.1 

4  7 
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•  r. 
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*♦  u 
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fa  1 

4T 
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H.«. 
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4  * 

6  1 

M  4 
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4 

4  *- 
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*♦  7 
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4  7 

fa  7 

4  H 

4 

4  r 

J  V 
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t  r, 

49 

5  “ 

'  0 
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r.  c 

fc  7 

- 1 

.1 

L  1 

:  <  2 

r  | 

>  c. 

6  7 

^  3 

.14 

.  _  X 

(i  ‘ 

..  4 

Hi 

v  l 

»-  , 

■  74 

j  7 

■■  b 

1  f- 

;  7 

6  i- 

7 

>7 

1  & 

S-) 

.fa 

8  i1' 

b  9 
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‘.9 

39 

i.O 

fal 

6  0 

7  0 

i,l 

6  2 

r  1 

71 

t  7 

hi 

>2 

72 

63 

7  7 

Lb 

fcn 
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69 
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71 
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7b 
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76 
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7  ft 
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79 
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60 

12« 

Hi 
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H3 
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b  4 
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H4 
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64 

13  3 

Hb 
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ri  b 
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B  fa 

1 3fa 

7  b 
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13fc 

76 
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£8 

1  4  C 

77 
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6  0 
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14  1 

-7  0 
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78 

14b 

“1 

2  4  fa 

70 
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53 
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94 

14  0 

94 

lbo 

•j  b 
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^0 
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9  h 
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nb 
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97 

1  bb 

9  D 

1  bt 

bh 
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9  7 
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5  c- 

lbC 
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>1 
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66 

79 

70 

8  0 

70 
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70 
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70 
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64 
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73 
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74 
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7b 

87 

76 
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76 

88 

77 
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77 

89 

66 

90 

68 
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7b 

91 

78 

91 

79 
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79 
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»0 

94 

a3 

84 

83 

9b 

84 

8  5 
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96 

bb 

66 

bb 

97 

8  fa 

9  8 
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9  e 

8  7 

99 
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68 
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88 
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49 
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b9 
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“3 
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91 
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1 
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1 
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1 
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COORDINATES  FOR  THE  70.00  F  ISOTHERM 


0.0227 

-0.2658 

0.0000 

-0.2705 

0.1032 

-0.2493 

0.1230 

-0.2355 

0.1892 

-0.1892 

0.2030 

-0.1668 

0.2436 

-0.1009 

0.2481 

-0.0715 

0.2886 

0.0000 

0.2802 

0.0331 

0.2334 

0.0967 

0.2101 

0.1279 

0.1746 

0.1746 

0.1378 

0.1971 

0.0928 

0.2241 

0.0481 

0.2324 

0.0000 

0.2410 

coordinates  for 

THE  40.00  F 

ISOTHERM 

0.0076 

-0.1366 

0.0536 

-0.1293 

0.0987 

-0.0987 

0.1284 

-0.0532 

0.1382 

0.0000 

0.126° 

0.0526 

u • C  96  6 

0 .0966 

0.0820 

0.1256 

0.0413 

U.1278 

O.COOC 

0.1358 

COOR  )  1 N A T  E  S  FOR 

THt  100.00  F 

ISOTHERM 

0.0383 

-0.0924 

0.0707 

“  L<  •  U  7  U  f 

0 . 0  °2  4 

-0.0383 

G.C924 

0  .  C  3  3 

0.0707 

0.0737 

0.0383 

3.0924 

G  •  o  0  C  D 

u  •  1  C  ’  0 
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Table  82.  Sample  input  to  FEHEAT-all  files  must  be  typed  in  by  user. 


a.  File  ED. 


L 

2 

3 

4 

5 
b 

7 

8 
9 

10 

11 

12 

13 

14 
lb 
lb 

17 

18 

19 

20 
21 
22 

23 

24 
2b 
26 

27 

28 

29 

30 

31 

32 

33 

34 
3b 
3b 

37 

38 

39 

40 

41 

42 
4  3 
44 
4b 

46 

47 

48 

49 

50 

51 
5? 
5* 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 
6b 
66 

67 

68 

69 

70 

71 

72 

73 

74 
7b 

76 

77 

78 

79 

80 
81 
82 
83 


10 

10 

11 

11 

12 

12 

13 

13 

14 

14 
13 

15 

1  o 

16 
17 
17 
19 
19 
23 
23 
?1 
21 
22 
22 
?  3 

23 

24 

24 

25 

25 

26 
26 
28 
28 

29 

30 

30 

31 

31 
3? 

32 

33 

33 
3  4 

34 

35 
35 

37 
3  7 

38 

38 

39 

3  9 

4  0 

4 : 

4  1 
41 
4  2 
4? 

43 
4  3 

44 
44 
46 

46 

47 

47 

48 

48 

49 
44 
3C 

50 

51 

51 

52 
5? 

53 
53 
55 

55 

56 


? 

11 

3 

12 

4 

13 
6 

14 
6 

15 

7 

16 

8 
17 

q 

Id 

11 

20 

12 

21 

13 
22 

14 

23 

15 

24 

16 

25 

17 

26 

18 
27 
20 

29 
21 

30 
22 

31 

23 

32 

24 
3  3 

25 
3  4 

26 
’  5 
2  7 
3 1 
2  9 
36 
30 


1 

2 

2 

3 

3 

4 

4 

5 

5 

6 
6 
7 

7 

8 
8 
9 

10 

11 

11 

12 

12 

13 

13 

14 

14 

15 

15 

16 
16 
17 


31 

4  3 

32 
4  1 

33 
42 

34 
4  3 

35 

44 

36 

45 
36 

47 
39 

48 

4  0 

49 

41 

50 

42 

5 1 

43 

52 

44 

5  3 

*5 

54 

47 
56 

48 


18 

19 

20 
29 
21 
?1 
22 
22 
23 

23 

24 

24 

25 

25 

26 
?  6 
?7 
26 
29 

29 
*0 

30 

31 

31 

52 

32 

33 

53 

34 

34 

35 
3b 
3b 

37 

38 

38 

39 
59 

40 

40 

41 

41 

42 

42 

43 

43 

44 

44 

45 
4b 
47 
47 


1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

l 

1 

1 

1 


1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


84 

85 
8b 

87 

88 

89 

90 

91 

92 

93 

94 

95 
9b 

97 

98 

99 
100 
101 
102 

103 

104 

105 
10b 

107 

108 

109 

110 
HI 
112 

113 

114 

115 

116 

117 

118 

119 

120 
121 
122 

123 

124 

125 

126 

127 

128 

129 

130 

131 

132 
135 

134 

135 
13b 
137 

13  6 

139 

140 

14  1 

142 

143 

144 

145 
14b 

147 

148 

149 

150 

151 

152 

153 

154 

155 

156 

157 

158 

159 

160 
161 
162 

163 

164 

165 

166 


5b 

57 

57 

58 

58 

59 

59 

60 
60 
61 
61 
62 
62 
64 

64 

65 
65 
b  & 
f->  £> 
67 

67 
bft 

68 

69 
s9 

70 

70 

71 
71 
64 
73 
66 
77 
bft 

73 
70 
32 
7  3 

74 

75 

76 
76 
79 
‘•O 
81 


ft  9 


e  4 


84 


8  5 


6  6 

7b 

8  7 

76 
3  8 

77 
8° 

6  b 

9  C 

78 
ft  1 
74 
ii 
a4 
o4 
95 

95 
9b 

96 

97 
86 

98 

87 

99 

88 
100 

89 
R9 

90 

90 

91 
91 


57 

49 

58 

50 

59 

51 

60 

52 
61 

53 
62 

54 
63 

56 
65 

57 
6b 

58 

67 

59 

68 

60 

69 
61 

70 
62 

71 
63 

72 

73 
6b 

77 
68 

78 
70 
82 
72 

74 

75 

76 

77 
7° 

8  0 
B  1 
8  2 

84 

73 

74 

85 

75 
ftb 
8  7 
87 
P.  ft 
b  ft 
89 
89 


5  0 

91 

91 

92 
')? 

94 
B  4 

95 

85 
9  6 

86 
97 
9ft 
9ft 
99 
99 

100 

100 

101 

101 

102 

102 

103 

103 

104 


48 

48 

49 

49 

50 
60 

51 

51 

52 

52 

53 

53 

54 

55 

56 

56 

57 

57 

58 

58 

59 

59 

60 
60 
61 
61 
62 
62 
63 

65 
b5 
67 
67 

69 
b9 
71 
71 

66 
66 
r  6 
6  b 

70 
70 
70 
70 
64 
64 

73 

74 
7« 

75 
75 
7b 
7b 
77 

77 
68 
48 
73 

78 

79 
79 
HC 
83 

83 

34 

84 

65 

85 

86 
86 
8  7 

87 

88 
88 
89 

89 
102 

90 

103 

91 

104 
«2 


1 

1 

l 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

l 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 
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Table  B2  (cont'd). 
b.  File  NPD. 


1 

C  .0  3  C 0 

-0.1003 

2 

3.0383 

-  0  . !!  9  2  4 

5 

3.0707 

-  e  .  0  7  2  7 

4 

l<  .  C  )24 

-2  .  1"  3 

5 

3.1000 

c.woeo 

b 

C  .  0  9  1 4 

3.3383 

7 

0.3707 

0.2707 

b 

C  •  C  3 1>  3 

0.0924 

9 

3 . 0  3  a  c 

C  .  1  0  3  0 

10 

C  .02  20 

-0.1503 

1 1 

0.0574 

-0. 1386 

1  2 

C.1361 

-  2 . 1  0  6 1 

13 

0 . 1  3u  6 

-2.2574 

19 

3.1500 

0  .  C  0  C  0 

15 

2.1366 

3.2574 

16 

C  .  1  0  6 1 

C.1C61 

1  7 

0.0674 

0.1366 

18 

0.03CC 

0.1500 

19 

0.0030 

-C.2C0C 

2  C 

0.3765 

-0.164  !- 

21 

'1.14  14 

-0.1414 

22 

3. 1&4& 

-  C . O  7 65 

23 

3.2200 

0.0000 

29 

0 . 1 1'  4  b 

0.0765 

2b 

0.1414 

0.1414 

26 

C .0765 

n  •  1 6  4  8 

27 

2.0000 

0.3000 

28 

2.0030 

-2.5033 

29 

0.114b 

-  0  •  '  7  7  2 

30 

3.2121 

—0.2121 

31 

3.2772 

-0.114'< 

32 

2.3030 

3.000C 

33 

0.2772 

0.1149 

39 

0.2121 

0.2121 

35 

3.1148 

0.2772 

36 

C.0000 

0.3030 

37 

3.0000 

-0.42 00 

38 

0.1637 

-0.1883 

39 

2.2970 

-t.i  °70 

9  C 

C  ,36.80 

—  C  .  1  6  2  7 

9  1 

2.4203 

0  •  0  0  0  3 

92 

3.1883 

C.16C7 

93 

3 . 2 Q  73 

C . 2  n  70 

44 

3.1607 

0.76  H r 

4  5 

0.030C 

0.4200 

46 

0.0033 

-0.38CC 

47 

3.2220 

-0.5359 

46 

[.MCI 

-0.4101 

49 

0.5359 

-0.2220 

50 

0.560C 

0.0020 

51 

0.5359 

0.2220 

52 

0.4101 

C  .  4  1  0  1 

53 

0.2220 

0.5359 

54 

j  .  0  o  C  0 

0.5602 

55 

o.coco 

-0.7802 

56 

2. 29' 5 

-0. 72C6 

57 

C  .5bl5 

-0.6515 

58 

6.7236 

-0.2985 

59 

3 .7300 
0.72 2 6 

C.OOCO 

60 

0.2965 

61 

0.5515 

0.5515 

62 

3.2985 

0.723b 

63 

C.C003 

0.7800 

64 

0.3320 

-1. 0003 

65 

3.3327 

-0.9239 

66 

C.7C71 

-0.7071 

67 

0.9235 

-3.3827 

ba 

l.orcc 

3.20PC 

69 

3.6239 

0.  3827 

70 

3.7C71 

0.7071 

71 

0 . IP  2  7 

0.«?S9 

72 

a .  o  o  e  o 

1.0030 

7  3 

3.4142 

-l.COOO 

74 

0.7071 

-1.0000 

75 

1.0GC3 

-1.0000 

76 

1.0300 

-a. 7071 

77 

1.0320 

-3.4142 

7b 

1.00C0 

0.4142 

79 

1.0000 

0.7071 

80 

1.0330 

l.OOCO 

51 

0.7371 

l.COOO 

32 

0.4142 

1.0000 

33 

C  .  0  3  0  0 

-1.4200 

84 

2.4300 

-1.4302 

35 

2.6330 

-1.43PC 

86 

1.4303 

-1.4000 

87 

1 .400  0 

-1.0000 

8  8 

1.4C00 

-C.fcOCO 

89 

1 . 4  0  C  0 

- C.20CD 

93 

1.4DC0 

0.2000 

91 

1 .42100 

0.6000 

92 

1.4000 

1.300C 

93 

0  .  C  0  3  0 

-2.3000 

94 

3.4020 

-2.1000 

95 

3.8220 

-2.0000 

96 

1.4020 

-2.30CO 

97 

2.0000 

-2.0000 

98 

2.C000 

-1.4000 

99 

2.0000 

-l.COOO 

10  0 

2.0000 

-C.600C 

131 

2.0000 

-0.20C0 

102 

2.0330 

0.2000 

103 

2.0002 

0.6000 

104 

2.000C 

l.OOCO 

82 


c.  File  BCT. 


1 

1 

10Q.0CQU 

2 

'l 

1  00 .00  OC 

A 

S 

ICO. 0000 

<* 

<* 

1  0  0 . 0  C  0  c 

5 

5 

100.0000 

6 

6 

10O.3CC0 

7 

7 

ICO.OCOO 

8 

8 

ICO.CjOO 

S 

9 

ICO. 00 C  3 

10 

HO 

10.0000 

1  1 

Hi 

10.0003 

12 

n? 

10. coco 

IS 

72 

10.0003 

M 

93 

10.0000 

15 

10.0000 

16 

9  5 

1C. 0300 

1  7 

h6 

10.0000 

18 

9  7 

10.0000 

19 

H8 

10.0000 

20 

99 

1C. 0000 

21 

100 

10.0300 

22 

1  01 

10.0000 

23 

1  0  2 

10.0300 

2M 

1  C3 

1C . 03  OC 

25 

1  3* 

10.0000 

26 

92 

10.0000 

d.  File  QUAN. 

1C*  lbe  2  6  0  0  U  fe 


e.  File  TIOT. 

10.00 
SO. 00 
50.00 
TO. 00 
90.00 
100.00 
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A  facsimile  catalog  card  in  Library  of  Congress  MARC 
format  is  reproduced  below. 
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